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

Abstract

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.

Background

The prevalence of renal cell carcinoma (RCC), a lethal urogenital cancer, ranks third after prostate and bladder cancers (1–3). In 2020, about 73,750 new RCC cases were diagnosed, with approximately 14,830 deaths in the USA (3). Nowadays, a range of treatments, such as surgery accompanied with or without postoperative adjuvant therapy, chemotherapy, immunotherapy, and target therapy, have been developed for RCC. Although these options have certain therapeutic effects, the overall prognosis of RCC patients remains dismal, especially in the late-stage RCC (4).

Over recent decades, the development of immunotherapy has revolutionized cancer treatment paradigms and has been recognized as a promising therapeutic frontier (5–7). For example, immune checkpoint blockade (ICB) is a new therapeutic strategy for several cancer types, such as breast cancer (8, 9), melanoma (10, 11), and lung cancer (12, 13). ICB has also evolved in RCC and showed certain practical application value through the years based on the phase III CheckMate 025 study, whether or not patients have been previously treated (14, 15). In addition, accumulating evidence has also proven that the tumor immune microenvironment (TIME), which encompasses immune cells, fibroblasts, extracellular matrix, endothelial cells, and various cytokines, is associated with tumor progression and metastasis (16–19). In 2017, Chevrier et al. depicted an in-depth Immune Atlas of Clear Cell Renal Cell Carcinoma by applying mass cytometry for the high-dimensional single-cell analysis of kidney primary Tumors (20). In addition, an increased number of studies have proved that multiple immune cells, including CD8⁺ T cells, CD4⁺ T cells and NK cells et al, have been associated with ccRCC tumor (21, 22). An in-depth understanding of TIME is critical to identifying potential immunotherapeutic targets for RCC. However, the majority of the studies have only evaluated gene expressions in the prediction of survival rates for RCC patients, and most of these biomarkers only reveal the status of TIME in some aspects (23, 24). Hence, a comprehensive immune-based model might provide an in-depth insight into the association between prognosis and TIME in RCC.

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.

Materials and Methods

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.

Establishment of the Immune-Related Risk Score

Among 538 KIRC with mRNA expression data, 517 patients with the overall survival (OS) data were retained for further analyses. First, 70% of samples were randomly drawn and grouped as training cohort to develop a prognostic risk model, and the other 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” R-function 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 immune-associated 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.

As chemotherapy and targeted therapy are widely used to treat clear cell renal cell carcinoma (ccRCC), risk scores were used to predict the drug sensitivity based on half-maximal inhibitory concentrations (IC50) for each KIRC patient from the Genomics of Drug Sensitivity in Cancer (GDSC) website (38) using the R package “pRRophetic” (39–44). The common target drugs, such as Cisplatin, Gefitinib, Gemcitabine, Sorafenib, Sunitinib, Vinblastine, Vinorelbine, and Vorinostat, were selected for ccRCC.

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. R 4.0.3 and SPSS 26.0 software were used for all analyses.

Results

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).

Figure 1

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).

Figure 2

Establishment and Validation of Prognostic Immune Score Model

All 367 KIRC samples in the training cohort were utilized in a prognostic model establishment. First, univariate Cox regression analysis was carried out to explore the association between DE-IRGs and the OS outcomes for KIRC samples. Among 98 DE-IRGs, 47 genes were selected. To avoid overfitting, we further conducted the LASSO Cox regression analysis with minimized lambda (Figures 3A, B). A total of 14/47 genes were used to establish the prognostic immune score model using the following formula: risk score = (SAA1 × 0.08215) + (IL20RB × 0.07643) + (TNFSF14 × 0.09743) + (ESRRG × -0.09743) + (FGF21 × 0.23324) + (IFNG × 0.05956) + (CTLA4 × 0.01439) + (KLRK1 × 0.00717) + (IL11 × 0.01639) + (GDF6 × - 0.02484) + (BMP7 × 0.05433) + (GNLY × 0.11780) + (AVPR1B × -0.06460) + (CXCL11 × 0.03907) (Figure 3C).

Figure 3

Based on the above calculated formula, the risk scores for every patient in the training set were computed. Then, with the median risk score as the basis, patients were allocated into the high- and low-risk groups. The high-risk patients exhibited significantly poor OS outcomes compared to low-risk patients (P=1.398E-10) (Figure 4A). As illustrated in Figure 4B, the AUCs of risk scores for 1-, 3-, and 5-years were 0.725, 0.723, and 0.745, respectively, in the training group. The distribution of the risk scores, survival time, survival status, and the expression of 14 OS-associated DE-IRGs for KIRC patients in the training cohort are displayed in Figures 4C–E. The model’s C-index was 0.698 (95% confidence interval (CI): 0.647–0.750, P=6.849E-14). To further explore whether the prognostic model was independent of other clinical elements, such as grade, age, T stage, and clinical stage, univariate and multivariate Cox regression analyses were conducted (Table 2). The risk score was confirmed as an independent prognostic factor (HR=2.699, 95% CI: 1.716–4.243, P<0.0001).

Figure 4

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 5C–E and Figures 6C–E. The C-indices of the model were 0.835 (95% CI: 0.782–0.888, P=1.034E-35) and 0.709 (95% CI: 0.666-0.752, P=7.520E-22) in the testing and entire cohorts, respectively. Univariate and multivariate Cox regression analyses for clinicopathological parameters were carried out in the testing 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). The C-indices of the nomogram in the testing and entire cohorts were 0.859 (95% CI: 0.809–0.909, P=4.980E-45) and 0.786 (95% CI: 0.752-0.821, P=4.201E-59), respectively. Finally, the calibration curves and decision curves for 1-, 3-, and 5-year survival probabilities were established (Figures 5H, I and Figures 6H, I). These findings indicated that the nomogram has excellent predictive performance in all cohorts.

Figure 5

Figure 6

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 Figures 7B, C, respectively. Compared to the low-risk group, the high-risk score group exhibited significantly elevated proportions of plasma cells, T cells CD8⁺, T cells follicular helper, T regulatory cells (Tregs), and M0 macrophages (P<0.05) (Figures 7C, D). Conversely, the proportions of macrophages M1, activated natural killer (NK) cells, naïve B cells, macrophages M2, resting NK cells, monocytes, T cells CD4⁺ memory resting, and resting mast cells in the high-risk score subtype were remarkably elevated compared to those in the low-risk score subtypes (P<0.05) (Figures 7C, D). In addition, in 22 immune cell types, high plasma cells, Tregs, follicular helper T cells, and monocytes M0 level were remarkably correlated with 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.

Figure 7

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–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.

Figure 8

To further characterize immune cell infiltration, 28 immune cell signatures (25, 49–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 immune-associated 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 tumor-specific CD⁺8 T cells by activating CXCL10 (54–57). Some studies (57–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_infiltration_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 TP53-related pathways exhibited significantly elevated NES values in the high-risk patient group, whereas the Hippo-, NRF2-, PI3K-, RAS-, and TGF-β-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-β-related pathways had a negative correlation with the risk score (Figure S8B).

Therapeutic Benefit of the Risk Score

Recently, ICB therapies have exhibited striking clinical benefits. However, the main challenge faced by ICB therapies is the limitation of effective predictive markers with only a few patients showing therapeutic response. Herein, the urothelial cancer database (IMvigor210) consisting of anti-PD-L1 therapy and the malignant melanoma database (GSE91016) administered with anti-PD-1 and-CTLA-4 therapy were used to investigate the association between risk score and immunotherapeutic benefits. Figures 9A–F and Figures S9A–C showed the distribution of clinical and molecular characteristics (immunotherapy response, binary response, immune phenotype, immune cells (IC) level, and tumor cells (TC) level and correlation with risk scores 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 chi-squared 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 low-risk 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).

Figure 9

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, CTLA4_pos_PD1_neg, CTLA4_neg_PD1_pos, and 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).

Figure 10

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).

Figure 11

Risk Score and TMB

Next, we analyzed the gene mutations of each KIRC patient. The waterfall chart showed the top 20 genes with the highest mutation frequencies: VHL, PBRM1, SETD2, MTOR, TTN, MUC16, KDM5C, BAP1, HMCN1, DNAH9, LRP2, ATM, ARID1A, CSMD3, DST, KMT2C, ERBB4, SMARCA4, USH2A, 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).

Figure 12

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.

Figure 13

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 π-π stacking interactions. Hence, these two compounds were both regarded as potential SAA1 inhibitors that could improve the prognosis of RCC.

Table 4

Name	Compound Structure	Target	Binding Energy (kcal/mol)	Combination Type
merbromin		SAA1	-7.85	Hydrogen bonds, Hydrophobic interactive, π-stacking
cefamandole		SAA1	-7.43	Hydrogen bonds, Hydrophobic interactive, π-stacking
fludrocortisone		SAA1	-7.35	Hydrogen bonds, Hydrophobic interactive
cetirizine		SAA1	-7.26	Hydrogen bonds, Hydrophobic interactive, π-stacking
finasteride		SAA1	-7.09	Hydrogen bonds, Hydrophobic interactive
vincamine		SAA1	-6.96	Hydrogen bonds, Hydrophobic interactive
sulfathiazole		SAA1	-6.01	Hydrogen bonds, Hydrophobic interactive
biperiden		SAA1	-5.61	Hydrophobic interactive, π-stacking
natamycin		SAA1	-5.35	Hydrophobic interactive, π-stacking
orphenadrine		SAA1	0	0

The selected compounds of docking results.

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 immune-escape 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–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–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 antigen-presenting cell (APC) recruitment, activation, and antigen-specific immune responses (78). CTLA-4 is a negative regulator and modulates T cell activation, and induces tolerance (79). CXCL11 is activated by IFN-γ and IFN-β 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–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.

Funding

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

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1

  SiegelRLMillerKDJemalA. Cancer Statistics, 2019. CA Cancer J Clin (2019) 69:7–34. doi: 10.3322/caac.21551

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  SiegelRLMillerKDJemalA. Cancer Statistics, 2020. CA Cancer J Clin (2020) 70:7–30. doi: 10.3322/caac.21590

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  SiegelRLMillerKDFuchsHEJemalA. Cancer Statistics, 2021. CA: A Cancer J Clin (2021) 71:7–33. doi: 10.3322/caac.21654

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  MakhovPJoshiSGhataliaPKutikovAUzzoRGKolenkoVM. Resistance to Systemic Therapies in Clear Cell Renal Cell Carcinoma: Mechanisms and Management Strategies. Mol Cancer Ther (2018) 17:1355–64. doi: 10.1158/1535-7163.MCT-17-1299

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  AounFRassyEEAssiTAlbisinniSKatanJ. Advances in Urothelial Bladder Cancer Immunotherapy, Dawn of a New Age of Treatment. Immunotherapy (2017) 9:451–60. doi: 10.2217/imt-2017-0007

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  HoosA. Development of Immuno-Oncology Drugs — From CTLA4 to PD1 to the Next Generations. Nat Rev Drug Discov (2016) 15:235–47. doi: 10.1038/nrd.2015.35

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  KamalYChengCFrostHRAmosCI. Predictors of Disease Aggressiveness Influence Outcome From Immunotherapy Treatment in Renal Clear Cell Carcinoma. Oncoimmunology (2019) 8:e1500106. doi: 10.1080/2162402X.2018.1500106

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BuXYaoYLiX. Immune Checkpoint Blockade in Breast Cancer Therapy. Adv Exp Med Biol (2017) 1026:383–402. doi: 10.1007/978-981-10-6020-5_18

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  HuZIHoAYMcArthurHL. Combined Radiation Therapy and Immune Checkpoint Blockade Therapy for Breast Cancer. Int J Radiat Oncol Biol Phys (2017) 99:153–64. doi: 10.1016/j.ijrobp.2017.05.029

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  MahoneyKMFreemanGJMcDermottDF. The Next Immune-Checkpoint Inhibitors: PD-1/PD-L1 Blockade in Melanoma. Clin Ther (2015) 37:764–82. doi: 10.1016/j.clinthera.2015.02.018

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  Van AllenDMiaoDSchillingBShuklaSABlankCZimmerLet al. Genomic Correlates of Response to CTLA-4 Blockade in Metastatic Melanoma. Science (2015) 350:207–11. doi: 10.1126/science.aad0095

  • CrossRef
  • Google Scholar

• 12

  AnagnostouVSmithKNFordePMNiknafsNBhattacharyaRWhiteJet al. Evolution of Neoantigen Landscape During Immune Checkpoint Blockade in Non-Small Cell Lung Cancer. Cancer Discov (2017) 7:264–76. doi: 10.1158/2159-8290.CD-16-0828

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  XiaLLiuYWangY. PD-1/PD-L1 Blockade Therapy in Advanced Non-Small-Cell Lung Cancer: Current Status and Future Directions. Oncologist (2019) 24:S31–41. doi: 10.1634/theoncologist.2019-IO-S1-s05

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  MotzerRJEscudierBMcDermottDFGeorgeSHammersHJSrinivasSet al. Nivolumab Versus Everolimus in Advanced Renal-Cell Carcinoma. N Engl J Med (2015) 373(19):1803–13. doi: 10.1056/NEJMoa1510665

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  MotzerRJTannirNMMcDermottDFArén FronteraOMelicharBChoueiriTKet al. Nivolumab Plus Ipilimumab Versus Sunitinib in Advanced Renal-Cell Carcinoma. N Engl J Med (2018) 378(14):1277–90. doi: 10.1056/NEJMoa1712126

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  CooperLAGutmanDAChisolmCAppinCKongJRongYet al. The Tumor Microenvironment Strongly Impacts Master Transcriptional Regulators and Gene Expression Class of Glioblastoma. Am J Pathol (2012) 180:2108–19. doi: 10.1016/j.ajpath.2012.01.040

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  CurryJMSprandioJCognettiDLuginbuhlABar-adVPribitkinEet al. Tumor Microenvironment in Head and Neck Squamous Cell Carcinoma. Semin Oncol (2014) 41:217–34. doi: 10.1053/j.seminoncol.2014.03.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  SenbabaogluYGejmanRSWinerAGLiuMVan AllenEMde VelascoGet al. Tumor Immune Microenvironment Characterization in Clear Cell Renal Cell Carcinoma Identifies Prognostic and Immunotherapeutically Relevant Messenger RNA Signatures. Genome Biol (2016) 17:231. doi: 10.1186/s13059-016-1092-z

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  YoshiharaKShahmoradgoliMMartinezEVegesnaRKimHTorres-GarciaWet al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  ChevrierSLevineJHZanotelliVRTSilinaKSchulzDBacacMet al. An Immune Atlas of Clear Cell Renal Cell Carcinoma. Cell (2017) 169(4):736–749.e18. doi: 10.1016/j.cell.2017.04.016

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  KopeckýOLukesováSVroblováVVokurkováDMorávekPSafránekHet al. Phenotype Analysis of Tumour-Infiltrating Lymphocytes and Lymphocytes in Peripheral Blood in Patients With Renal Carcinoma. Acta Med (Hradec Kralove) (2007) 50(3):207–12. doi: 10.14712/18059694.2017.84

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  KomoharaYHasitaHOhnishiKFujiwaraYSuzuSEtoMet al. Macrophage Infiltration and Its Prognostic Relevance in Clear Cell Renal Cell Carcinoma. Cancer Sci (2011) 102(7):1424–31. doi: 10.1111/j.1349-7006.2011.01945.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  WangBLiuBYuGHuangYLvC. Differentially Expressed Autophagy-Related Genes Are Potential Prognostic and Diagnostic Biomarkers in Clear-Cell Renal Cell Carcinoma. Aging (2019) 11(20):9025–42. doi: 10.18632/aging.102368

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  WangYChenLWangGChengSQianKLiuXet al. Fifteen Hub Genes Associated With Progression and Prognosis of Clear Cell Renal Cell Carcinoma Identified by Coexpression Analysis. J Cell Physiol (2019) 234:10225–37. doi: 10.1002/jcp.27692

  • CrossRef
  • Google Scholar

• 25

  BhattacharyaSAndorfSGomesLDunnPSchaeferHPontiusJet al. ImmPort: Disseminating Data to the Public for the Future of Immunology. Immunol Res (2014) 58:234–9. doi: 10.1007/s12026-014-8516-1

  • CrossRef
  • Google Scholar

• 26

  MariathasanSTurleySJNicklesDCastiglioniAYuenKWangYet al. Tgfβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature (2018) 554:544–8. doi: 10.1038/nature25501

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  LoveMIHuberWAndersS. Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2. Genome Biol (2014) 15:550. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  YuGWangLGHanYHeQY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  GroeneveldCSChagasVSJonesSJMRobertsonAGPonderBAJMeyerKBet al. RTNsurvival: An R/Bioconductor Package for Regulatory Network Survival Analysis. Bioinformatics (2019) 35:4488–9. doi: 10.1093/bioinformatics/btz229

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  RT. The Lasso Method for Variable Selection in the Cox Model. Stat Med (1997) 16:385–95. doi: 10.1002/(SICI)1097-0258(19970228)16:4<385::AID-SIM380>3.0.CO;2-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  IasonosASchragDRajGVPanageasKS. How to Build and Interpret a Nomogram for Cancer Prognosis. J Clin Oncol (2008) 26:1364–70. doi: 10.1200/JCO.2007.12.9791

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  DuanJXieYQuLWangLZhouSWangYet al. A Nomogram-Based Immunoprofile Predicts Overall Survival for Previously Untreated Patients With Esophageal Squamous Cell Carcinoma After Esophagectomy. J Immunother Cancer (2018) 6:100. doi: 10.1186/s40425-018-0418-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  KiranMChatrathATangXKeenanDMDuttaA. A Prognostic Signature for Lower Grade Gliomas Based on Expression of Long Non-Coding RNAs. Mol Neurobiol (2019) 56:4786–98. doi: 10.1007/s12035-018-1416-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  VickersAJElkinEB. Decision Curve Analysis: A Novel Method for Evaluating Prediction Models. Med Decis Making (2006) 26:565–74. doi: 10.1177/0272989X06295361

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  NewmanAMLiuCLGreenMRGentlesAJFengWXuYet al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  RooneyMSShuklaSAWuCJGetzGHacohenN. Molecular and Genetic Properties of Tumors Associated With Local Immune Cytolytic Activity. Cell (2015) 160:48–61. doi: 10.1016/j.cell.2014.12.033

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  CharoentongPFinotelloFAngelovaMMayerCEfremovaMRiederDet al. Pan-Cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  Geeleher PCNHuangRS. Clinical Drug Response can be Predicted Using Baseline Gene Expression Levels and In Vitro Drug Sensitivity in Cell Lines. Genome Biol (2014) 15(3):R47. doi: 10.1186/gb-2014-15-3-r47

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  LiuZZhangYShiCZhouXXuKJiaoDet al. A Novel Immune Classification Reveals Distinct Immune Escape Mechanism and Genomic Alterations: Implications for Immunotherapy in Hepatocellular Carcinoma. J Transl Med (2021) 19(1):5. doi: 10.1186/s12967-020-02697-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  LiuZWangLGuoCLiuLJiaoDSunZet al. TTN/OBSCN ’Double-Hit’ Predicts Favourable Prognosis, ’Immune-Hot’ Subtype and Potentially Better Immunotherapeutic Efficacy in Colorectal Cancer. J Cell Mol Med (2021) 25(7):3239–51. doi: 10.1111/jcmm.16393

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  LiuZLiuLLuTWangLLiZJiaoDet al. Hypoxia Molecular Characterization in Hepatocellular Carcinoma Identifies One Risk Signature and Two Nomograms for Clinical Management. J Oncol (2021) 2021:6664386. doi: 10.1155/2021/6664386

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  LiuZLiuLJiaoDGuoCWangLLiZet al. Association of RYR2 Mutation With Tumor Mutation Burden, Prognosis, and Antitumor Immunity in Patients With Esophageal Adenocarcinoma. Front Genet (2021) 12:669694. doi: 10.3389/fgene.2021.669694

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  LiuZWangLLiuLLuTJiaoDSunYet al. The Identification and Validation of Two Heterogenous Subtypes and a Risk Signature Based on Ferroptosis in Hepatocellular Carcinoma. Front Oncol (2021) 11:619242. doi: 10.3389/fonc.2021.619242

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  LiuZLuTWangLLiuLLiLHanX. Comprehensive Molecular Analyses of a Novel Mutational Signature Classification System With Regard to Prognosis, Genomic Alterations, and Immune Landscape in Glioma. Front Mol Biosci (2021) 8:682084. doi: 10.3389/fmolb.2021.682084

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  MayakondaALinDCAssenovYPlassCKoefflerHP. Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res (2018) 28:1747–56. doi: 10.1101/gr.239244.118

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  AnXZhuYZhengTWangGZhangMLiJet al. An Analysis of the Expression and Association With Immune Cell Infiltration of the cGAS/STING Pathway in Pan-Cancer. Mol Ther - Nucleic Acids (2019) 14:80–9. doi: 10.1016/j.omtn.2018.11.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  LiTChengHYuanHXuQShuCZhangYet al. Antitumor Activity of cGAMP via Stimulation of cGAS-cGAMP-STING-IRF3 Mediated Innate Immune Response. Sci Rep (2016) 6:19049. doi: 10.1038/srep19049

  • CrossRef
  • Google Scholar

• 48

  RamanjuluJMPesiridisGSYangJConchaNSinghausRZhangS-Yet al. Design of Amidobenzimidazole STING Receptor Agonists With Systemic Activity. Nature (2018) 564:439–43. doi: 10.1038/s41586-018-0705-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  BechtEGiraldoNALacroixLButtardBElarouciNPetitprezFet al. Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression. Genome Biol (2016) 17:218. doi: 10.1186/s13059-016-1070-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  BindeaGMlecnikBTosoliniMKirilovskyAWaldnerMObenauf AnnaCet al. Spatiotemporal Dynamics of Intratumoral Immune Cells Reveal the Immune Landscape in Human Cancer. Immunity (2013) 39:782–95. doi: 10.1016/j.immuni.2013.10.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  MiaoYRZhangQLeiQLuoMXieGYWangHet al. ImmuCellAI: A Unique Method for Comprehensive T-Cell Subsets Abundance Prediction and Its Application in Cancer Immunotherapy. Adv Sci (Weinh) (2020) 7:1902880. doi: 10.1002/advs.201902880

  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  NirmalAJReganTShihBBHumeDASimsAHFreemanTC. Immune Cell Gene Signatures for Profiling the Microenvironment of Solid Tumors. Cancer Immunol Res (2018) 6:1388–400. doi: 10.1158/2326-6066.CIR-18-0342

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  WolfDMLenburgMEYauCBoudreauAvan ’t VeerLJ. Gene Co-Expression Modules as Clinically Relevant Hallmarks of Breast Cancer Diversity. PloS One (2014) 9:e88309. doi: 10.1371/journal.pone.0088309

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  Hildner KEBPurthaWE. Batf3 Deficiency Reveals a Critical Role for CD8alpha+ Dendritic Cells in Cytotoxic T Cell Immunity. Science (2008) 322:1097. doi: 10.1126/science.1164206

  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  FuertesMBKachaAKKlineJWooS-RKranzDMMurphyKMet al. Host Type I IFN Signals Are Required for Antitumor CD8+ T Cell Responses Through CD8α+ Dendritic Cells. J Exp Med (2011) 208:2005–16. doi: 10.1084/jem.20101159

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  RobertsEWBrozMLBinnewiesMHeadleyMBNelsonAEWolfDMet al. Critical Role for CD103(+)/CD141(+) Dendritic Cells Bearing CCR7 for Tumor Antigen Trafficking and Priming of T Cell Immunity in Melanoma. Cancer Cell (2016) 30:324–36. doi: 10.1016/j.ccell.2016.06.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  SprangerSDaiDHortonBGajewskiTF. Tumor-Residing Batf3 Dendritic Cells Are Required for Effector T Cell Trafficking and Adoptive T Cell Therapy. Cancer Cell (2017) 31:711–23.e4. doi: 10.1016/j.ccell.2017.04.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 58

  BöttcherJPBonavitaEChakravartyPBleesHCabeza-CabrerizoMSammicheliSet al. NK Cells Stimulate Recruitment of Cdc1 Into the Tumor Microenvironment Promoting Cancer Immune Control. Cell (2018) 172:1022–37.e14. doi: 10.1016/j.cell.2018.01.004

  • CrossRef
  • Google Scholar
  • View reference in article

• 59

  SprangerSBaoRGajewskiTF. Melanoma-Intrinsic β-Catenin Signalling Prevents Anti-Tumour Immunity. Nature (2015) 523:231–5. doi: 10.1038/nature14404

  • CrossRef
  • Google Scholar
  • View reference in article

• 60

  HarlinHMengYPetersonACZhaYTretiakovaMSlingluffCet al. Chemokine Expression in Melanoma Metastases Associated With CD8+ T-Cell Recruitment. Cancer Res (2009) 69:3077–85. doi: 10.1158/0008-5472.CAN-08-2281

  • CrossRef
  • Google Scholar
  • View reference in article

• 61

  ChengWCTsuiYCRagusaSKoelzerVHMinaMFrancoFet al. Uncoupling Protein 2 Reprograms the Tumor Microenvironment to Support the Anti-Tumor Immune Cycle. Nat Immunol (2019) 20:206–17. doi: 10.1038/s41590-018-0290-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 62

  Ogunleye AAWQ. XGBoost Model for Chronic Kidney Disease Diagnosis. IEEE/ACM Trans Comput Biol Bioinform (2019) 17(6):2131–40. doi: 10.1109/TCBB.2019.2911071

  • CrossRef
  • Google Scholar
  • View reference in article

• 63

  ZhongJSunYPengWXieMYangJTangX. XGBFEMF: An XGBoost-Based Framework for Essential Protein Prediction. IEEE Trans Nanobioscience (2018) 17:243–50. doi: 10.1109/TNB.2018.2842219

  • CrossRef
  • Google Scholar
  • View reference in article

• 64

  FerlayJColombetMSoerjomataramIDybaTRandiGBettioMet al. Cancer Incidence and Mortality Patterns in Europe: Estimates for 40 Countries and 25 Major Cancers in 2018. Eur J Cancer (2018) 103:356–87. doi: 10.1016/j.ejca.2018.07.005

  • CrossRef
  • Google Scholar
  • View reference in article

• 65

  FerlayJColombetMSoerjomataramIMathersCParkinDMPinerosMet al. Estimating the Global Cancer Incidence and Mortality in 2018: GLOBOCAN Sources and Methods. Int J Cancer (2019) 144:1941–53. doi: 10.1002/ijc.31937

  • CrossRef
  • Google Scholar
  • View reference in article

• 66

  MotzerRJJonaschEMichaelsonMDNandagopalLGoreJLGeorgeSet al. NCCN Guidelines Insights: Kidney Cancer, Version 2.2020. J Natl Compr Cancer Network (2019) 17:1278–85. doi: 10.6004/jnccn.2019.0054

  • CrossRef
  • Google Scholar
  • View reference in article

• 67

  BerglundAAmankwahEKKimYCSpiessPESextonWJManleyBet al. Influence of Gene Expression on Survival of Clear Cell Renal Cell Carcinoma. Cancer Med (2020) 9:8662–75. doi: 10.1002/cam4.3475

  • CrossRef
  • Google Scholar
  • View reference in article

• 68

  KnottMEMinattaJNRouletLGueglioGPasikLRanuncoloSMet al. Circulating Fibroblast Growth Factor 21 (Fgf21) as Diagnostic and Prognostic Biomarker in Renal Cancer. J Mol biomark Diagn (2016) 1(Suppl 2):015. doi: 10.1158/1538-7445.AM2016-431

  • CrossRef
  • Google Scholar
  • View reference in article

• 69

  Lin JYMXuXWangYXingHAnJYangJet al. Identification of Biomarkers Related to CD8 T Cell Infiltration With Gene Co-Expression Network in Clear Cell Renal Cell Carcinoma. Aging (Albany NY) (2020) 12(4):3694–712. doi: 10.18632/aging.102841

  • CrossRef
  • Google Scholar
  • View reference in article

• 70

  MarkicDCelicTGrskovicASpanjolJFuckarZGrahovacBet al. mRNA Expression of Bone Morphogenetic Proteins and Their Receptors in Human Renal Cell Carcinoma. Urol Int (2011) 87:353–8. doi: 10.1159/000330797

  • CrossRef
  • Google Scholar
  • View reference in article

• 71

  PanDXuLLiuHZhangWLiuWLiuYet al. High Expression of Interleukin-11 Is an Independent Indicator of Poor Prognosis in Clear-Cell Renal Cell Carcinoma. Cancer Sci (2015) 106:592–7. doi: 10.1111/cas.12638

  • CrossRef
  • Google Scholar
  • View reference in article

• 72

  XuFGuanYZhangPXueLYangXGaoKet al. The Impact of TNFSF14 on Prognosis and Immune Microenvironment in Clear Cell Renal Cell Carcinoma. Genes Genomics (2020) 42:1055–66. doi: 10.1007/s13258-020-00974-0

  • CrossRef
  • Google Scholar

• 73

  ParetCSchönZSzponarAKovacsG. Inflammatory Protein Serum Amyloid A1 Marks a Subset of Conventional Renal Cell Carcinomas With Fatal Outcome. Eur Urol (2010) 57(5):859–66. doi: 10.1016/j.eururo.2009.08.014

  • CrossRef
  • Google Scholar
  • View reference in article

• 74

  CuiXFCuiXGLengN. Overexpression of Interleukin-20 Receptor Subunit Beta (IL20RB) Correlates With Cell Proliferation, Invasion and Migration Enhancement and Poor Prognosis in Papillary Renal Cell Carcinoma. J Toxicol Pathol (2019) 32:245–51. doi: 10.1293/tox.2019-0017

  • CrossRef
  • Google Scholar
  • View reference in article

• 75

  KangMHChoiHOshimaMCheongJHKimSLeeJHet al. Estrogen-Related Receptor Gamma Functions as a Tumor Suppressor in Gastric Cancer. Nat Commun (2018) 9:1920. doi: 10.1038/s41467-018-04244-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 76

  VenkatesanAMVyasRGramannAKDresserKGujjaSBhatnagarSet al. Ligand-Activated BMP Signaling Inhibits Cell Differentiation and Death to Promote Melanoma. J Clin Invest (2017) 128:294–308. doi: 10.1172/JCI92513

  • CrossRef
  • Google Scholar
  • View reference in article

• 77

  RauletDH. Roles of the NKG2D Immunoreceptor and Its Ligands. Nat Rev Immunol (2003) 3:781–90. doi: 10.1038/nri1199

  • CrossRef
  • Google Scholar
  • View reference in article

• 78

  TewaryPYangDde la RosaGLiYFinnMWKrenskyAMet al. Granulysin Activates Antigen-Presenting Cells Through TLR4 and Acts as an Immune Alarmin. Blood (2010) 116:3465–74. doi: 10.1182/blood-2010-03-273953

  • CrossRef
  • Google Scholar
  • View reference in article

• 79

  ChambersCAKMichaelSEgen, JacksonGAllison, JamesP. CTLA-4-Mediated Inhibition in Regulation of T Cell Responses Mechanisms and Manipulation in Tumor Immunotherapy. Annu Rev Immunol (2001) 19(1):565–94. doi: 10.1146/annurev.immunol.19.1.565

  • CrossRef
  • Google Scholar
  • View reference in article

• 80

  TokunagaRZhangWNaseemMPucciniABergerMDSoniSet al. CXCL9, CXCL10, CXCL11/CXCR3 Axis for Immune Activation – A Target for Novel Cancer Therapy. Cancer Treat Rev (2018) 63:40–7. doi: 10.1016/j.ctrv.2017.11.007

  • CrossRef
  • Google Scholar
  • View reference in article

• 81

  FuchsCSDoiTJangRWMuroKSatohTMachadoMet al. Safety and Efficacy of Pembrolizumab Monotherapy in Patients With Previously Treated Advanced Gastric and Gastroesophageal Junction Cancer: Phase 2 Clinical KEYNOTE-059 Trial. JAMA Oncol (2018) 4:e180013. doi: 10.1001/jamaoncol.2018.0013

  • CrossRef
  • Google Scholar
  • View reference in article

• 82

  PandaAMehnertJMHirshfieldKMRiedlingerGDamareSSaundersTet al. Immune Activation and Benefit From Avelumab in EBV-Positive Gastric Cancer. J Natl Cancer Inst (2018) 110:316–20. doi: 10.1093/jnci/djx213

  • CrossRef
  • Google Scholar
  • View reference in article

• 83

  RohWChenPLReubenASpencerCNPrietoPAMillerJPet al. Integrated Molecular Analysis of Tumor Biopsies on Sequential CTLA-4 and PD-1 Blockade Reveals Markers of Response and Resistance. Sci Transl Med (2017) 9(379):eaah3560. doi: 10.1126/scitranslmed.aah3560

  • CrossRef
  • Google Scholar
  • View reference in article

• 84

  GoodmanAMKatoSBazhenovaLPatelSPFramptonGMMillerVet al. Tumor Mutational Burden as an Independent Predictor of Response to Immunotherapy in Diverse Cancers. Mol Cancer Ther (2017) 16:2598–608. doi: 10.1158/1535-7163.MCT-17-0386

  • CrossRef
  • Google Scholar
  • View reference in article

• 85

  KorbakisDSoosaipillaiADiamandisEP. Study of Kallikrein-Related Peptidase 6 (KLK6) and Its Complex With Alpha1-Antitrypsin in Biological Fluids. Clin Chem Lab Med (2017) 55:1385–96. doi: 10.1515/cclm-2017-0017

  • CrossRef
  • Google Scholar
  • View reference in article