Immunogenomic Analyses of the Prognostic Predictive Model for Patients With Renal Cancer

Background Renal cell carcinoma (RCC) is associated with poor prognostic outcomes. The current stratifying system does not predict prognostic outcomes and therapeutic benefits precisely for RCC patients. Here, we aim to construct an immune prognostic predictive model to assist clinician to predict RCC prognosis. Methods Herein, an immune prognostic signature was developed, and its predictive ability was confirmed in the kidney renal clear cell carcinoma (KIRC) cohorts based on The Cancer Genome Atlas (TCGA) dataset. Several immunogenomic analyses were conducted to investigate the correlations between immune risk scores and immune cell infiltrations, immune checkpoints, cancer genotypes, tumor mutational burden, and responses to chemotherapy and immunotherapy. Results The immune prognostic signature contained 14 immune-associated genes and was found to be an independent prognostic factor for KIRC. Furthermore, the immune risk score was established as a novel marker for predicting the overall survival outcomes for RCC. The risk score was correlated with some significant immunophenotypic factors, including T cell infiltration, antitumor immunity, antitumor response, oncogenic pathways, and immunotherapeutic and chemotherapeutic response. Conclusions The immune prognostic, predictive model can be effectively and efficiently used in the prediction of survival outcomes and immunotherapeutic responses of RCC patients.

In the current study, we established an immune prognostic signature model for RCC using the training cohort and further confirmed the effectiveness of the prognosis model in the testing and the entire cohort. Additionally, the associations between the risk score subtypes and immune checkpoints, antitumor immunity, antitumor response, oncogenic pathways, immune cell infiltration, and tumor mutation burden (TMB) were explored. Also, the models' ability for the prediction of chemotherapeutic and immunotherapeutic responses was evaluated. Finally, we screened out two compounds that could improve the prognosis of RCC.

Data Acquisition as Well as Preprocessing
Transcriptional expression profiles, mutation patterns, and related clinical data for KIRC patients were retrieved from the Cancer Genome Atlas (TCGA) cohort (https://cancergenome. nih.gov/). Immune-associated genes (IRGs) were derived from the Immunology Database as well as Analysis Portal (ImmPort) database (25). The immunophenoscore (IPS) for RCC patients were retrieved from the cancer immune group atlas (TCIA) (https://tcia.at/home). In addition, the advanced urothelial cancer database of administered anti-PD-L1 therapy was downloaded using the R package "IMvigor210CoreBiologies" (version 1.0.0) (26). The malignant melanoma dataset that received anti-PD-1 and antiCTLA4 therapy were obtained from the GSE91061 cohort. All data were subjected to background correction and logarithmic conversion using R software.

Differentially Expressed Immune-Related Genes (DE-IRGs) and Functional Enrichment Analyses
Differential gene expression analysis between tumor and corresponding normal tissues in KIRC were screened based on the count data for TCGA kidney cancer cohort using the R package "DESeq2" (27), according to the screening criteria (log2| fold change| >2, P-value <0.05). The IRGs involved in oncogenesis were provide by IMMPORT website. Then, DE-IRGs were identified by the intersection between DEGs and IRGs.
The R package "clusterProfiler" was used for Gene Ontology (GO) as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of these significant DE-IRGs and their visualization (28). Next, we defined the pathways and terms using false discovery rate (FDR) ≤0.05 as statistically significant.
30% of samples comprised the validation set, which was used in evaluating the model's predictive ability and robustness in the entire cohort. Then, DE-IRGs were screened out by univariate Cox proportional hazard regression through the "coxph" Rfunction from the "survival" package (29). Subsequently, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was carried out to select the prognostic genes using the R package "glmnet" (30). Finally, the immuneassociated risk score was calculated using LASSO Cox regression hazard regression − retrieved regression coefficients to multiply expression levels of genes (the risk score = mRNA expression levels of gene a × coefficient a + mRNA expression levels of gene b × coefficient b + ……+ mRNA expression levels of gene n × coefficient n).
In addition, by setting the median of risk score as the cutoff value, the patients were classified into a high-risk group and a low-risk group. To establish the prognostic accuracy of the established model, we used Kaplan-Meier survival curve analysis, concordance (C)-index, log-rank test in addition to time-dependent receiver operating characteristic curves (ROC) and XGBoost algorithm.

Independent Prognostic Value of the Immune-Associated Prognostic Signature
Multivariate Cox regression analysis with the forward stepwise procedure was performed to investigate if the risk score is an independent prognostic factor. The immune-associated risk score and other clinical variables with P <0.05 were identified as independent prognostic risk factors.

Establishment and Validation of the Nomogram
To develop a prognostic signature for 1-, 3-, and 5-year survival rates, a nomogram was constructed using the identified independent prognostic variables, such as stage, age, and risk score (31). Moreover, the C-index, calibration curve, decision curve analysis (DCA), and ROC analysis were performed to determine its predictive accuracy and discriminatory capacity (32). The C-index was evaluated using a bootstrap method involving 1000 resamples (33). The C-index values, dependent on the nomogram's predictive ability, ranged from 0.5 (no discrimination) to 1 (perfect discrimination). The consistency between the predictive survival rate and the actual survival rate in unknown samples was assessed using calibration curves. Additionally, DCA (34) was used to evaluate the clinical utility and the net benefits of the nomogram as it takes both discrimination and calibration into consideration. Finally, the area under the receiver operating characteristic (ROC) curve (AUC) was also determined for each variable to evaluate the discriminative performance of the nomograms.

Immune Cell Proportion Analyses and Immune Related Features
To explore immune cell abundance in KIRC tissues, CIBERSORT (35) was employed to evaluate the proportions of 22 immune cell types using a deconvolution algorithm by the R package with default parameters. In addition, the ESTIMATE scores (ES), tumor purity (TP), stromal scores (SS), and immune scores (IS) for each KIRC sample were evaluated using the ESTIMATE algorithm (19) of the "estimate" package. The cytolytic activity (CYT) index is a geometric mean of mRNA expression levels of GZMA and PRF1, and was utilized to assess the intratumoral immune cytolytic T-cell activities (36).
Immunotherapy and Chemotherapeutic Response in Risk Score Subtype As immune checkpoint molecules are widely explored in the immunotherapeutic studies of multiple cancers, programmed cell death 1 (PDCD1, also referred to as PD-1), CD274 molecule (also referred to as PD-L1), and cytotoxic T-lymphocyte protein 4 (CTLA4) were used to evaluate the associations between risk scores and immunotherapeutic efficacies. The urothelial cancer dataset (IMvigor210) comprising of administered anti-PD-L1 therapy was used to establish the therapeutic benefits between high-and low-risk score subtypes using four treatment categories: progressive disease (PD), stable disease (SD), complete response (CR), and partial response (PR).
IPS is a machine learning-based scoring system applied for the prediction of patients' responses to immune checkpoint inhibitor (ICI) treatment based on the weight average Z scores representing immune-related genes expression in cell types (37). High IPS scores reflect increased immunogenicity.
Tumor Mutational Burden (TMB), Connectivity Map (CMAP) and Molecular Docking Analysis KIRC patients' somatic variants data were analyzed and visualized by "maftool" R package (45) to identify the mutation burden of KIRC in the high-and low-risk scores. Then, the TMB of each patient was calculated as follows: mutations/million bases.
Next, to identify the potentially small molecules related to this signature, genes in the model were assessed via CMAP analysis. Thus, the positive mean represented that these selected drugs may share similar functions with the model, while the negative mean indicated that these drugs could improve the prognosis of RCC. Herein, we screened compounds by setting the criteria as P <0.05.
Moreover, the crystal structure of the protein was obtained from RCSB Protein Data Bank. The three-dimensional structures for all compounds were downloaded from PubChem database using MOL2 format. The molecular docking calculations were conducted using Schrodinger and Pymol 2.1 software.

Statistical Analysis
The differences between variables were determined by chi-square as well as Student's t-tests. For baseline clinical data, the Wilcoxon test and the Kruskal-Wallis were utilized to evaluate the significant differences between two or multiple groups, respectively. The Kaplan-Meier survival curves were compared using the log-rank test. P<0.05 indicated statistical significance.

Identification and Functional Analyses of DE-IRGs
All 517 KIRC samples with OS information were split into training (367 patients) and test groups (151 patients). Between the training and validation cohorts, no significant differences were detected among most of the clinical characteristics ( Table 1).
With the cutoff value |log2 fold change (logFC)|>2 and adjusted P<0.05, 953 DEGs were filtered, of which 539 genes were significantly elevated, while 414 genes were significantly suppressed in tumor samples compared to normal samples ( Figures 1A, B). Moreover, principal component analysis (PCA) results ( Figure 1C) revealed that KIRC samples clustered separately from normal samples. Subsequently, the intersection between DEGs and immune-associated genes retrieved from the ImmPort database was determined, and 98 DE-IRGs were selected and visualized on a Venn diagram ( Figure 1D).
These 98 DE-IRGs were further utilized in functional enrichment analyses, including KEGG and GO analyses. Based on GO analysis, in the biological process (BP), these DE-IRGs were enriched in cell chemotaxis, leukocyte chemotaxis, lymphocyte chemotaxis, positive regulation of cell adhesion, and T cell activation ( Figure 2A). In the cellular component (CC) category, the DE-IRGs were mainly enriched in the cytoplasmic vesicle lumen, plasma membrane's external side, platelet alpha granule, platelet alpha granule lumen, and vesicle lumen ( Figure 2B). Regarding molecular function (MF), these DE-IRGs were enriched in cytokine receptor binding, growth factor activity, receptor ligand activity, cytokine activity, and signaling receptor activator activity ( Figure 2C). Regarding KEGG pathways analysis, these DE-IRGs were mainly involved in the calcium signaling pathway, chemokine signaling pathways, cytokine-cytokine receptor interactions, Ras signaling pathways, and viral protein interactions with cytokine receptors and cytokines ( Figure 2D).
In addition, a quantitative strategy for the prediction of the prognostic outcomes of patients was established by constructing a nomogram that integrated the risk scores as well as other independent clinical prognostic factors for OS ( Figure 4F). Then, the nomogram's performance was determined using the ROC curve, C-index, calibration curve, and decision curve analyses. The AUCs of the nomogram were 0.828, 0.783, and 0.774 for 1-, 3-, and 5-year survival times, respectively ( Figure 4G). The C-index was 0.762 (95% CI: 0.720-0.804, P=1.800E-34). Based on the calibration curve, the training cohort predicted that 1-, 3-, and 5-year survival probabilities were good ( Figure 4H). For the decision curve, the nomogram exhibited a higher net benefit than other schemes to predict the OS ( Figure 4I).
To delineate the robustness and versatility of the immune score model, the risk score in the training cohort was validated in the testing and entire cohorts. The participants in the testing and entire cohorts were grouped into high-and low-risk score subtypes using the same formula. The findings in the testing and entire datasets were similar. The Kaplan-Meier survival curves revealed poor survival rates for the high-risk group in the testing (P=3.58E-8) ( Figure 5A) and the entire cohorts (P=3.616E-12) ( Figure 6A). The AUC for 1-, 3-, and 5-years are 0.858, 0.842, and 0.857 in the testing group ( Figure 5B) and 0.736, 0.727, and 0.746 in the entire group ( Figure 6B). The survival data, risk score, scatterplots, and gene expression pattern distributions in the testing and entire cohorts are shown in Figure  and entire cohorts. Also, the risk score was an independent prognostic indicator of OS in KIRC patients ( Table 2). To improve the prognostic immune score model, the nomogram system was established based on testing and entire cohorts ( Figure 5F and Figure 6F). The AUC of our nomogram for predicting 1-, 3-, and 5-year OS was 0.9, 0.875, and 0.891, respectively, in the testing cohort and 0.858, 0.808, and 0.787, respectively, in the entire cohort ( Figure 5G and Figure 6G).

Associations Between DE-IRGs Signature and Clinical Characteristics of KIRC Patients
Next, we further investigate the association between clinical characteristics, including tumor burden, age at diagnosis, gender, grade, clinical stage, T stage, and the prognostic risk signature. A significant correlation was established between high-risk score and a high tumor burden (P=9.87E-08), male gender (P=0.03), advanced grade (P=2.54E-09), higher stage (P=8.85E-11), and T stage (P=5.6E-08) ( Figure S1). Additionally, no statistical significance was observed between < 65-year-old group and >65-year-old group (P=0.1). Subsequently, we also assessed whether the model could assess the survival probability in subgroups exhibiting varying clinical patterns. The prognostic model could be utilized for the prediction of survival probabilities for various clinicopathological parameters (P<0.05) ( Figure S2).

Immune Cell Proportions Between High-and Low-Risk Score Patients
Using the CIBERSORT algorithm, 22 immune cell types were determined in each KIRC sample between high-and low-risk score subtypes. The proportions of 22 immune cells and their distribution in tumor samples are illustrated in Figure 7A and poor OS outcomes (P=0.01, 0.0019, <0.0001, and 0.031, respectively), while the increase in activated dendritic cells was related to better OS (P=0.0079) ( Figure S3). Figure S4 displayed a weak or moderate correlation between the levels of various tumor-infiltrating immune cells and the risk score.

Immune Landscape in KIRC Patients
Subsequently, the associations between risk score and some immune-associated features were assessed. The cGAS-STING pathway has been shown to be a key signaling pathway in antitumor immunity and cancer therapeutics (46)(47)(48). Thus, four key genes (TBK1, IRF3, MB21D1, and TMEM173) in the cGAS-STING signaling pathway, three immune checkpoint molecules (PD-L1, CTLA-4, and PD-1), CYT, and the results of ESTIMATE algorithm (SS, IS, ES, and TP) and risk score were investigated. Figure 8A shows that the risk score values are correlated with the immune score, tumor purity, TBK1, ESTIMATE score, IRF3, stromal score, MB21D1, PD-1, and CTLA-4. Figure S5 showed significant differences in the CYT, immune score, ESTIMATE score, stromal score, and tumor purity based on the Wilcoxon test between the two risk score subtypes (P<0.0001). Importantly, the expression of IRF3, MB21D1, TMEM173, PD-1, and CTLA-4 was elevated in the high-risk than in the low-risk score subtype.
To further characterize immune cell infiltration, 28 immune cell signatures (25,(49)(50)(51)(52)(53) from diverse resources were investigated based on the single sample gene set enrichment analysis (ssGSEA) algorithm. As shown in Figure 8B, 23 immune subpopulations (multiple T cell signatures, including T helper cells, central memory CD8+ T cells, and activated CD T cells) were enriched in high-risk patient cohort, whereas only two subpopulations (immature dendritic cells and neutrophils) were enriched in the low-risk patient group. Furthermore, DEGs between low-and high-risk groups were determined by gene set enrichment analysis (GSEA) using two MeSH terms (gene2pubmed and gendoo) to explore their immune-related functions. The DEGs were enriched in multiple immuneassociated terms, including CD4-CD8 ratio, immune tolerance, lymphocyte cooperation, lymphocyte count, immunologic memory, and T-cell antigen receptor specificity in gendoo and gene2pubmed ( Figures 8C, D).

Correlation Between Risk Score Model and T Cell Infiltrations, Antitumor Immunity, Antitumor Responses, and Oncogenic Pathways
Several studies have shown that cDC1 cells play a central role in the initiation of antitumor CD8 + T cells and driving tumorspecific CD + 8 T cells by activating CXCL10 (54)(55)(56)(57). Some studies (57)(58)(59) also clarified that the two key chemokines (CCL4 and CCL5) are the key modulators of cDC1 recruitment into tumors via activating CCR5 expression. Moreover, chemokines CXCR3, CXCL9, and CXCL10 have been documented on T cell infiltration and NK cell recruitment (60). Thus, we investigated the expression level of CCL4, CXCR3, CXCL9, CCL5, and CXCL10 between high-and low-risk subtypes and the correlations between these genes and the risk score. The high-risk group patients exhibited higher expression levels compared to low-risk patients (P<0.05) ( Figures S6A-E). Moreover, strong positive correlations were established between risk scores and CXCR3, CCL5, CXCL9, CCL4, and CXCL10 (P<0.05) (Figures S6F-J).
Moreover, we explored the association between risk scores, T cell infiltrations, and antitumor response scores (BATF3_DC, IFNA, IFNG, IL_1_speed, T cell_infi ltration_1, T cell_infiltration_2, TFNA, and TNFa_speed) determined by ssGSEA from the corresponding TME gene signatures (57,61). For the high-risk group, the ssGSEA scores for T cell infiltrations and antitumor responses were significantly elevated compared to the low-risk group, as determined by the Wilcoxon test (P<0.05) ( Figure S7A). A strong positive correlation was established between risk scores and ssGSEA scores save to BATF3_DC (P<0.05) (Figures S7B-I). Conclusively, high-risk score patients exhibited elevated T cell infiltration levels.
The differences in the normalized enrichment score (NES) value of 10 oncogenic pathways between low-and high-risk groups were calculated using ssGSEA algorithm; also, the correlation between the NES value and the risk score was evaluated. Compared to the low-risk group, cell cycle and TP53related pathways exhibited significantly elevated NES values in the high-risk patient group, whereas the Hippo-, NRF2-, PI3K-, RAS-, and TGF-b-related pathways in the high-risk patient group had lower NES value (P<0.05) ( Figure S8A). The correlations between the risk score and the NES value in the cell cycle (P=1.44e-12) and TP53-related (P=0.024) pathways were found to be positive ( Figure S8B). Nevertheless, we also observed that the NES value of the Hippo-, NRF2-, PI3K-, RAS-, and TGF-b-related pathways had a negative correlation with the risk score ( Figure S8B).  between high-and low-risk groups in the IMvigor210 cohort and GSE91061 cohort separately. For the immunotherapy response, the risk score of RCC with CR/PR were significantly lower than those of RCC with SD/PD, as assessed by the chisquared test (IMvigor210 dataset: P<0.001, GSE91061 dataset: P=0.036) ( Figure 9A and Figure S9C). The violin plot further revealed that the risk scores in CR/PR were lower than those in SD/PD, as assessed by the Wilcoxon test (IMvigor210 cohort: P=1.3e-08, GSE91061 cohort: P=0.0075) ( Figure 9C and Figure  S9B). Strikingly, Kaplan-Meier curves showed that high-risk score patients exhibited worse prognosis compared to the lowrisk score patients in IMvigor210 (P<0.0001) ( Figure 9G) and GSE91061 cohort (P=0.00016) ( Figure S8D). In addition, IMvigor210 and GSE91061 were used to plot a time-dependent ROC. The current results displayed that the AUCs of our model for OS were 0.61 at 6 months, 0.673 at 12 months, and 0.729 at 18 months in the IMvigor210 cohort ( Figure 9H) and 0.746 at 12 months, 0.712 at 18 months, and 0.753 at 24 months in the GSE91061 cohort ( Figure S9A).
To further expand this study, the machine learning-based score (IPS) was determined to predict patients' response to ICI treatment. Four subtypes of IPS values (CTLA4_neg_PD1_neg, C T L A 4 _ p o s _ P D 1 _ n e g , C T L A 4 _ n e g _ P D 1 _ p o s , an d CTLA4_pos_PD1_pos) were carried out to predict the KIRC patients' responses to anti-CTLA4 and anti-PD1 treatment. We found that relative probabilities to response to anti-PD1 were elevated in high-risk score patients (P=0.023), and the similar results were obvious in the combination treatment of anti-PD1 and anti-CTLA4 (P=2.24e-04) ( Figure 10A). In addition, CTLA- 4 and PD-1 mRNA expression levels in the high-risk score group were significantly elevated compared to the low-risk score patients (P=1.07e-14 and P=2.02e-15), whereas no obvious difference was detected in the PD-L1 mRNA expression level between high-and low-risk patients (P=0.603) ( Figure 10B). This phenomenon was consistent with the concept that high expression of ICI genes had a poor prognosis. Owing to the complex environment between immune infiltration and ICI genes, we further examined whether immune infiltration had consequences on the clinical prognosis in ICI genes. Figure 10C shows that low-risk score patients with high PD-1 exhibited better clinical outcomes compared to high-risk score and high PD-1, and the outcomes of low-risk score patients with low PD-1 were superior to those of high-risk score patients and low PD-1 levels (P<0.0001). Also, patient groups showed similar findings, and survival patterns were yielded using risk score and PD-L1 or CTLA4 (P<0.0001) ( Figures 10C).
The responsive predictive values of the risk score to chemotherapy and target-therapy were also investigated by the IC50 of eight drugs. The estimated IC50 values of Cisplatin, Gemcitabine, Sorafenib, and Vinorelbine in high-risk patients were significantly elevated compared to low-risk patients, which indicating the high-risk patients showed a stronger drug resistance (P<0.05) (Figures 11A, B). Similarity, patients with high-risk group were associated with increased sensitivity to Gefitinib, Vinblastine, and Sunitinib relative to low-risk patients (P<0.05) (Figures 11A, B). and PCLO ( Figure 12A). Subsequently, the TMB for each sample was determined and was found to be higher in the high-risk patients (P=0.037) ( Figure 12B) and related to shorter OS (P=0.023) than in low-risk patients ( Figure 12C).

Prediction of High-and Low-Risk Scores by XGBoost Algorithm
XGBoost is an efficient and reliable machine learning classifier based on gradient boosting, designed to solve data science challenges accurately and rapidly in bioinformatics (62,63). Using this approach, a classifier that could predict high-and low-risk score groups for KIRC patients based on expression levels of 14 selected genes was constructed for the training cohort. SHAP dependency plot and the importance of 14 features were visualized in Figures 13A, B to evaluate the contribution of each feature towards prediction. Figure 13C showed that the AUC of the training cohort was 100%. Then, classification model performance was assessed using the testing and entire total cohorts ( Figures S10A, B and Figures 13C).
Taken together, the middle cutoff value might be suitable to classify KIRC patients.

Identification of Potential Small Molecule Drugs
According to CMAP analysis, 10 small molecule drugs with highly significant correlations are listed in Table 3. Among these, Finasteride, Biperiden, Merbromin, Cefamandole, Fludrocortisone, and Vincamine displayed a high negative correlation and potential to improve the prognosis of RCC. Subsequently, the SAA1 gene contributing to the model according to the feature importance was docked with these 10 compounds (Table 4). Next, we identified the compounds except for Orphenadrine that showed a high binding affinity against the target protein due to their binding energy <-5 kcal/mol. Moreover, the three-dimensional structure of top two high-affinity compounds combined with SAA1 is shown in Figures S11A, B.
In SAA1-merbromin complex, due to multiple phenylene rings and active groups, merbromin forms hydrogen bonds with activity groups of amino acids, such as GLN-66, ARG-25, and TRP-53, indicating that merbromin could match well with SAA1 protein.
Similarly, the SAA1-Cefamandole complex can be formed by multiple interactions, such as the cooperation of hydrogen bonding and multiple p-p stacking interactions. Hence, these two compounds were both regarded as potential SAA1 inhibitors that could improve the prognosis of RCC.

DISCUSSION
Epidemiological evidence indicated that the incidence of RCC had a continually increasing trend with high mortality (64,65). Clinical decision-making tools were effective prognostic biomarkers to predict the survival outcomes of RCC patients, rendering them a viable choice for clinicians. To date, the prognostic prediction of RCC patients relies on the TNM staging system according to the clinical practice guidelines (66). However, this system failed to taken the influence of gene level of RCC into consideration and made it not always able to predict the patients accurately. In recent years, IRGs have gradually gained attention with in-depth studies on immuneescape and immunotherapeutic mechanisms. Hence, an immune-related prognostic system is an urgent requirement for a supplementary TNM staging system. Next, we screened for immune-associated DEGs in RCC. To minimize the potential for overfitting, 14 genes established the prognostic immune signature and were validated in TCGA through the univariate Cox proportional hazard regression and LASSO Cox analysis. Subsequently, we confirmed the independent predictive roles of this signature. Then, a personalized, predictive nomogram with a risk score was developed, which served as a predictive indicator; the signature encompassed a total of 14 IRGs. Among these, SAA1, TNFSF14, FGF21, IFNG, BMP7, and IL11 are biomarkers for predicting RCC outcomes (67)(68)(69)(70)(71)(72). For example, as a member of the serum amyloid A family of apolipoproteins, SAA1 can increase the invasive capacity of tumor cells in RCC by inducing MMP-9 expression (73), which make it serve as a biomarker for the diagnosis and prognosis of advanced and metastatic renal cell carcinoma. In addition, as a member of the IL-6 family of cytokines, IL-11 exerts pleiotropic oncogenic activities may by stimulating angiogenesis and metastasis, which make it become an independent indicator of poor prognosis in RCC (71). The other IRGs, such as IL20RB, ESRRG, GDF6, were reported to be involved in the regulation of carcinogenesis (74)(75)(76) but not yet investigated in RCC. Moreover, some IRGs were also involved in TIME. For example, NKG2D receptor, KLRK1, is expressed in NK cells and activated CD8 + T cells, involved in innate immune responses (77). In some studies also identified GNLY as the first lymphocyte-derived alarmin protein to promote antigenpresenting cell (APC) recruitment, activation, and antigenspecific immune responses (78). CTLA-4 is a negative regulator and modulates T cell activation, and induces tolerance (79). CXCL11 is activated by IFN-g and IFN-b and can stimulate immune cells by promoting Th1 polarization and enhancing the antitumor immunity (80). To sum up, these IRGs may affected the prognosis and treatment of RCC by influencing TIME.
Herein, some self-validation processes, including the associations between risk scores and immune cell proportions, T cell infiltrations, antitumor immunity, antitumor response, GSEA analysis, and oncogenic pathways, were conducted to identify the risk score effectiveness in characterizing the immune landscape features of RCC patients. For immunotherapeutic development, anti-PD-1, anti-CTLA-4, and anti-PD-L1 treatment have been under intensive focus in solid tumors. Nevertheless, a small number of patients respond to such treatment, and some studies (81)(82)(83) pointed out that PD-L1 and PD-1 expression levels are not reliable biomarkers to predict ICI treatment. Hence, it is necessary for clinicians to develop a reliable tool for appropriate patient selection in immunotherapy. Based on these findings, we established that the risk score is a robust immune classifier for classifying RCC patients in different subtypes. Moreover, we also demonstrated that high-score patients were more immunotherapeutically suitable compared to patients in the low-risk score group. Targeted therapy is currently the main treatment strategy for metastatic RCC. Thus, it is necessary to identify patients with the potential to benefit from targeted therapy for RCC. Interestingly, our data showed that high-risk patients had a high sensitivity to Gefitinib, Vinblastine, and Sunitinib compared to low-risk score patients, who exhibited high sensitivity to Cisplatin, Sorafenib, Gemcitabine, and Vinorelbine. These responses could be attributed to the differences in the drug target. In addition, the TMB values of the high-risk score patients were elevated compared to those of the low-risk score patients. This finding was consistent with the concept that elevated TMB values are associated with a high probability of satisfactory immunotherapeutic outcomes (84,85).
Nevertheless, the present study had some limitations. First, although our model exhibited precise predictive capability to predict the survival of RCC patients, multiple large external cohorts of patients with RCC are also needed to further validate. Secondly, only the median risk score was used to classify the RCC patients into high-and low-risk score subtypes. An optimal cutoff of the risk score is essential for the stratification of RCC patients. Although our model had been correlated with immune cells, the mechanism underlying poor prognosis is unclear, requiring additional experimental and theoretical studies on immune cells in RCC to further understand their functional role.

CONCLUSIONS
Taken together, our proposed immune prognostic, predictive model could be used as a robust classifier for the prediction of survival outcomes and individual treatment guidance of adjuvant chemotherapy and anticancer immunotherapy for RCC.

DATA AVAILABILITY STATEMENT
The RNA-seq data and corresponding clinical information were observed from the TCGA (https://portal.gdc.cancer.gov/). The immune-related gene list was got from the IMMPORT website (https://www.immport.org/).

AUTHOR CONTRIBUTIONS
All authors participated in the design, interpretation of the studies, analysis of the data, and review of the manuscript. TF and JZ conceived and designed the whole project and wrote the manuscript. TF, DW, ZF, and ZW analyzed and visualized the data. TF, QL, PG, and XY interpreted the data and partook in the discussion. ML, YJ, and YL revised the final version of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the National Science Foundation of Beijing (7172068 and 7192053).