A model based on immunogenic cell death-related genes predicts prognosis and response to immunotherapy in kidney renal clear cell carcinoma

Introduction

Renal cell carcinoma (RCC) is one of the most common malignancies of the urinary system. RCC accounts for 2–3% of adult malignancies, and its incidence continues to increase (1). Currently, the incidence and mortality of RCC patients, particularly those with kidney renal clear cell carcinoma (KIRC), which represents 75% of all kidney cancers, has steadily improved (2). As the most common histological subtype of RCC, KIRC is known for its significant tumor heterogeneity, different clinical courses, and uncertainty in specific treatments (2,3).

With advances in medicine and surgery, cure rates and survival rates have improved greatly. Surgical resection remains the most common method for treating KIRC and cannot be replaced (4). By combining radiotherapy, chemotherapy, targeted therapy, and immunotherapy, many patients have been able to achieve longer survival rates (5,6). However, the molecular characteristics of KIRC have not yet been closely examined. To improve patient prognosis, the construction of reliable prognostic models might help guide clinical decisions in the treatment of KIRC patients. Previously, Chen et al. used three immune relative genes (IRGs) to develop a Cox regression model to predict the prognosis of patients with KIRC, and found that the model could accurately stratify patients with different survival outcomes (7). Wan et al. constructed a prognostic risk model (RM) using seven IRGs in renal papillary cell carcinoma and found that the model could independently distinguish between patients with different risks of death (8). However, few studies examining immunogenic cell death (ICD) in KIRC have been performed.

ICD can modulate cell death and activate the adaptive immune response of immunocompetent hosts (9). ICD includes the release of damage-associated molecular patterns (DAMPs), such as type-I interferons (IFNs), extracellular adenosine triphosphate (ATP), high-mobility group box-1 (HMGB1), heat shock proteins (e.g., HSP70 and HSP90), and the cell surface exposure of calreticulin (CRT). All ICDs originate from dying tumor cells and lead to the activation of tumor-specific immune responses (10,11). Using the immune system to trigger an anti-tumor immune response is a cancer immunotherapy concept. Numerous studies have shown the ability of ICDs to trigger certain anti-cancer immune response (12).

ICDs can activate tumor-specific immune responses through certain chemotherapeutic agents, oncolytic viruses, photodynamic therapy, and radiation therapy that stimulate the long-term therapeutic effects of anti-tumor drugs through a combination of direct cancer cell killing and anti-tumor immunity (12). For example, previous research has shown that radiotherapy and a few chemotherapeutic agents (including oxaliplatin and adriamycin) can cause the induction of ICD in vitro and in vivo and can also stimulate a tumor cell-specific immune response (13,14).

In light of the critical role that ICD plays in tumor progression, recent research has identified new ICD-related genetic signatures for the diagnosis or prognosis of cutaneous melanoma, head and neck squamous cell carcinoma, low-grade glioma, and breast cancer (11,15-17). In recent years, various pre-clinical studies have explored the molecular mechanism of ICD, but few studies have examined patients in a clinical context to evaluate the meanings of ICD (e.g., few studies have sought to identify biomarkers to classify patients according to their responses to ICD immunotherapy) (18). However, such research would be extremely beneficial. Thus, we first conducted a bioinformatics analysis to compare the expression levels of ICD-related genes between KIRC and adjacent normal tissues. We then assessed the prognostic value of the genes and analyzed the association between ICD and immune responses in the tumor microenvironment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-214/rc).

Methods

Data source

On October 1, 2022, the clinical parameters and RNA-sequencing (RNA-seq) data of 541 KIRC tissues and 72 adjacent normal tissues were retrieved from The Cancer Genome Atlas (TCGA) database. For the evaluation, we used the data provided by the International Cancer Genome Consortium (ICGC) project, specifically the Renal Cell Cancer-European Union (RECA-EU) data set, comprising 91 KIRC samples (https://dcc.icgc.org/projects/RECA-EU) with available RNA-seq data and clinical data. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

The differential expression profiles of the ICD-related genes

A total of 34 ICD-related genes were extracted from previous studies (19). Using the “limma” package, we identified the ICD-related genes that were differently expressed in the KIRC and the adjacent normal tissues (20). The identified protein names were uploaded to the functional protein association network, and protein-protein interaction (PPI) mapping was performed.

Formulation and verification of prognostic features

Univariate and multivariate Cox regression analyses were conducted to identify the prognosis-associated genes. Subsequently, we calculated the risk score using the following equation: risk score = Coef1 × Exp1 (Gene1) + Coef2 × Exp2 (Gene2) +... CoefN × ExpN (GeneN), where “coef” represents coefficient, and “exp” represents gene expression level.

Using the calculated median risk score as the cut-off value, 541 patients with KIRC from TCGA data set were allocated to the high- and low-risk groups. Next, using a Kaplan-Meier (K-M) analysis, we built the overall survival (OS) curves for the two groups. A receiver operating characteristic (ROC) analysis dependent on time was also conducted to examine survival. The predictive power of the model was assessed based on the area under the curve (AUC) of the ROC curve. Using risk scores and clinical data (e.g., age, sex, grade, and stage), nomogram models were constructed to predict the OS of KIRC patients. RECA-EU data were downloaded from the ICGC database for verification, and the survival differences between the different risk groups and the constructed prognosis model were verified.

Independent prognostic analysis of risk score

We downloaded patients’ clinical data (e.g., age, sex, grade, and stage) from TCGA data set, and we then constructed a univariate and multivariable Cox regression model to analyze the independent prognostic characteristics and risk scores of these variables.

Functional enrichment analysis

Initially, 956 differentially expressed genes (DEGs) between the high- and low-risk groups were identified in accordance with the following criteria: |log 2 fold change| >1.5 and false discovery rate <0.05). Subsequently, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Geno Ontology (GO) analyses were conducted.

Immune infiltration, gene variation, and drug-sensitivity analysis

We used the Estimation of STromal and Immune cells in MAlignant Tumors using Expression (ESTIMATE) algorithm to estimate the stromal and immune cells in the malignant tumor tissues and analyze the immune component and overall stromal component in the high- and low-risk subgroups of TCGA data set (21). The CIBERSORT algorithm was employed to evaluate the abundance of 22 immune cell types (22). A correlation analysis was also conducted to investigate the association between prognostic ICD-related genes and the immune infiltration landscape. Subsequently, the Wilcoxon test was used to compare the proportions of immune cells between the two groups. The Genomics of Drug Sensitivity in Cancer (GDSC) database was used to analyze the drug sensitivity of the ICD-related risk genes (23). Somatic mutation data of the KIRC samples in the “maf” format were obtained from TCGA GDC information portal. Next, using the “Maftools” (a package of R) in R software, waterfall plots were created to facilitate the visualization and aggregation of the mutated genes.

Statistical analysis

R software (version 4.1.3) was used to perform all the statistical analyses, and a P value <0.05 was considered statistically significant. The rank correlations among the different variables were assessed using the Pearson correlation coefficient test. Differences between the variables were assessed using the independent t-test.

K-M curves and log-rank tests were used to analyze the survival data, and a univariate Cox regression analysis was used to identify the factors affecting the survival of patients diagnosed with Clear cell renal cell carcinoma (ccRCC). A multivariate Cox regression analysis was used to identify the independent prognostic factors. A time-dependent ROC analysis was used to evaluate the accuracy of the prognostic prediction model. The model can be predicted according to this standard: low accuracy: 0.5< AUC-ROC ≤0.7, moderate accuracy: 0.7< AUC-ROC ≤0.9, and high accuracy: 0.9< AUC-ROC ≤1 (24).

Results

Identification of 2 ICD-related subtypes using a consensus clustering technique

Garg et al. conducted a comprehensive meta-analysis and summarized 34 ICD-related genes (19). The STRING database was used to analyze the PPI network to ascertain the correlations among the discovered ICD-related genes (Figure 1A). Subsequently, to better elucidate the relationships between these ICD-related genes, we built a co-expression network and investigated the associated functions (Figure 1B). ICD gene expression differed widely between the normal and KIRC samples (Figure 1C). Next, we used consensus clustering to identify two clusters in TCGA data set (Figure 1D-1F). A comparison of the majority of the expression levels of the ICD-related gene between clusters 1 and 2 revealed that the genes in cluster 1 were more highly expressed than those in cluster 2 (Figure 1G). Thus, we defined clusters 1 and 2 as the high- and low-ICD subtypes, respectively. The survival analysis revealed that the low-ICD subtype had a better prognosis than the high-ICD subtype (Figure 1H).

Figure 1 Determination of ICD-related subtypes through consensus clustering. (A) PPIs between ICD-related genes using the STRING database. (B) Analysis of DEGs and their co-expressed genes via GeneMANIA. (C) Heatmap of consensus clustering solutions (k =2) for 34 genes in 541 KIRC samples. (D-F) The delta area curves for the consensus clustering represent the comparative change in the area under the curve for the CDF curve for k =2 to 7. (G) Heatmap of the expression of the 34 ICD-related genes in distinct subtypes. Blue signifies lower expression, while red signifies higher expression. (H) Kaplan-Meier curves of overall survival in the cluster 1 and 2 subtypes. CDF, cumulative distribution function; ICD, immunogenic cell death; PPIs, protein-protein interactions; DEGs, differentially expressed genes; KIRC, kidney renal clear cell carcinoma.

RM construction of ICD-related genes

The univariate Cox regression analysis revealed that 10 ICD-related genes were linked to OS (Figure 2A). The ICD-related genes were further examined based on multivariate Cox regression analysis to build a RM. The formula for the risk score was expressed as follows: risk score = (0.17350 × LY96 exp.) + (0.11091 × FOXP3 exp.) + (−0.30199 × TLR4 exp.) + (0.27993 × IFNB1 exp.) + (0.09888 × IL6 exp.) (Figure 2B). After computing the risk scores for the KIRC patients, the median risk score was used as the threshold. The patients were then divided into high- and low-risk groups (Figure 2C).

Figure 2 Development and evaluation of prognostic features using ICD-related genes in TCGA data set. (A,B) Univariate Cox regression analysis to assess the prognostic significance of the ICD-related genes for overall survival. Distribution of risk scores from TCGA database. *, P<0.05; **, P<0.01; ***, P<0.001. (C) Each patient’s survival status and (D) Kaplan-Meier curve plots for high- and low-risk subgroups. (E,F) Heatmap of the prognostic five-gene profile. CI, confidence interval; ICD, immunogenic cell death; TCGA, The Cancer Genome Atlas.

We also focused on the relationship between survival statuses and risk scores. According to our findings, the low-risk group had a significant survival rate than the high-risk group (Figure 2D). Further, a K-M analysis was conducted to evaluate the prognostic significance of risk status in KIRC (Figure 2E). High-risk scores were found to be correlated with poor OS in TCGA cohort. A heat map of the results was generated and revealed distinctions in the expression of the five ICD-related genes between the high- and low-risk groups (Figure 2F). In the low-risk group, TLR4 was highly expressed, while LY96, FOXP3, IFNB1, and IL6 were lowly expressed. The survival rates of the ICD-related genes in the high- and low-expression groups were also used to construct K-M curves. Notably, low survival rates were observed, when LY96, FOXP3, IFNB1, and IL6 were highly expressed, or TLR4 was lowly expressed (Figure S1).

ICD risk score could be an independent predictor of OS

We conducted univariate and multivariate Cox regression analyses to demonstrate that the risk score obtained from the five-gene signature model might be an independent factor in prognosis. Both analyses confirmed that risk score was an independent prognostic factor that predicted poor survival in KIRC patients [univariate analysis hazard ratio (HR) =1.629 and multivariate analysis HR =1.626] (Figure 3A,3B). In addition, the ROC curves were plotted to validate the prognostic significance of the risk score (Figure 3C-3E). The AUC values were 0.725 for a 1-year period, 0.735 for a 3-year period, and 0.73 for a 5-year period, which shows that the risk score variable had a higher predictive value than the clinical variables. The RECA-EU data were downloaded from the ICGC database to verify these results, and the differences between the different risk groups and the constructed prognosis model was verified (Figure S2A-S2C). Additionally, using gene expression and the corresponding KIRC clinical data, we analyzed the correlations between the clinical variables and risk scores of the five ICD-related genes (i.e., LY96, FOXP3, IFNB1, IL6, and TLR4). Risk scores and FOXP3 were correlated with tumor grade, Metastasized stage (M stage), Regional Lymph Nodes stage (N stage), Tumor stage (T stage), and comprehensive stage. IFNB1 was associated with grade, M stage, T stage, and comprehensive stage; IL6 was associated with grade, T stage, and comprehensive stage. LY96 was associated with T stage and comprehensive stage; while TLR4 was associated with comprehensive stage (Figure S3A-S3T). In addition, we built a nomogram using the risk scores and clinical features to more precisely predict the prognosis of KIRC patients (Figure 4A,4B). Based on the ROC curve results, the nomogram was more accurate at predicting patient survival than the risk score alone (Figure 4C-4E). The results of both the univariate and multivariate Cox regression analyses showed that the nomogram was also an independent predictor of OS in KIRC patients (Figure 4F,4G).

Figure 3 Independent prognostic analysis of risk scores. (A) Prognostic impact analysis of KIRC risk scores and clinical characteristics using a univariate Cox regression analysis. (B) Independent prognostic impact analysis of risk scores and conventional prognostic clinical characteristics in KIRC using a multivariate Cox regression analysis. (C-E) Receiver operating characteristic analysis of KIRC risk scores and other prognostic clinical characteristics to predict the 1-, 3-, and 5-year survival rate of KIRC patients. CI, confidence interval; AUC, area under the curve; KIRC, kidney renal clear cell carcinoma.

Figure 4 Predictive nomogram creation. (A) Nomogram predicting the 1-, 2-, 3-, 4- and 5-year survival rates for KIRC patients. (B) Calibration curves for the nomogram prediction validation cohort of KIRC patient survival (x-axis: the predicted survival probability; y-axis: the actually observed survival probability). (C-E) Receiver operating characteristic analysis of the KIRC nomogram and other prognostic clinical characteristics to predict the 1-, 3-, and 5-year survival rates of KIRC patients. (F) Assessing the relationship between risk scores and clinical factors with OS using a forest plot of the univariate Cox test results. (G) The independent risk factors for OS were determined using a forest plot of the multivariate Cox analysis results. OS, overall survival; AUC, area under the curve; CI, confidence interval; KIRC, kidney renal clear cell carcinoma.

Identification of DEGs and signaling pathways in the low- and high-risk groups

To analyze the distinct mechanisms that influence the survival status of patients, we first identified 956 genes in TCGA data set that were differentially expressed between the two groups. We then conducted GO and KEGG functional enrichment analyses (Figure 5A,5B). We discovered that immunity and immune responses were differentially correlated with the low- and high-risk subgroups. Notably, according to the GO enrichment analysis results, the DEGs were enriched in immune-associated mechanisms, including the humoral immune response, immunoglobulin complex, and antigen binding (Figure 5C). The KEGG enrichment analysis demonstrated that the DEGs were significantly enriched in cytokine-cytokine receptor interactions, complement and coagulation cascades, and the IL-17 signaling pathways (Figure 5D).

Figure 5 Analysis of DEGs and analysis of enriched functions in distinct ICD risk scores. (A) Volcano map of DEG distribution in TCGA cohort’s ICD high- and low-risk groups. Red indicates that the ICD-related gene expression in the high-risk groups is higher than that in the low-risk groups. Green indicates that the ICD-related gene expression in the high-risk groups is lower than that in the low-risk groups. Black indicates that the difference between ICD-related gene expression in the high- and low-risk groups is not significant. (B) ICD high- and low-risk groups’ DEG expression in a heat map. (C,D) GO and KEGG enrichment analyses of signaling pathways; the different colors represent the q values. BP, biological process; CC, cellular component; MF, molecular function; DEGs, differentially expressed genes; ICD, immunogenic cell death; TCGA, The Cancer Genome Atlas; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Somatic mutations in the high- and low-risk groups

We first analyzed the genetic mutations in the high- and low-risk groups. The top five mutated genes in the low-risk group were von Hippel-Lindau (VHL) (46%), Polybromo 1 (PBRM1) (42%), Titin (TTN) (13%), SET domain containing 2 (SETD2) (9%), and mucin 16 (MUC16) (7%) (Figure 6A). The top five mutated genes in the high-risk group were VHL (46%), PBRM1 (39%), TTN (18%), SETD2 (16%), and BAP1 (15%) (Figure 6B). Compared to the low-risk group, the high-risk group had a higher incidence of TTN, SETD2, and BAP1 mutation rates.

Figure 6 The contrast of somatic mutations in distinct risk groups. (A,B) Visualization of the top 10 most commonly mutated genes in the high- (A) and low-risk groups (B) using Oncoprint.

Tumor microenvironment landscape in the high- and low-risk groups

Various immune estimation algorithms were used to further investigate the immune infiltration of patients with KIRC in the high- and low-risk groups. The findings of the ESTIMATE algorithm demonstrated that the high-risk group had decreased tumor purity, and improved immune scores, ESTIMATE scores, and stromal scores than the low-risk group (Figure 7A-7D). Subsequently, we used the “CIBERSORT” algorithm to assess the relative proportion of immune infiltration in 22 immune cell types and summarized the results of the KIRC patients from TCGA data sets (Figure 7E). The CIBERSORT algorithm results showed that the low-risk group had higher proportions of cluster of differentiation CD4⁺ resting memory T cells, resting natural killer (NK) cells, monocytes, M1 macrophages, M2 macrophages, resting dendritic cells, resting mast cells, activated mast cells, and eosinophils. Patients in the high-risk group had higher proportions of plasma cells, CD8⁺ T cells, activated memory CD4⁺ T cells, follicular helper T cells, regulatory T cells (Tregs), and M0 macrophages (Figure 7F). The single-sample gene-set enrichment analysis (ssGSEA) algorithm results revealed that 24 immune cells in the high-risk group scored significantly lower than those in the low-risk group, demonstrating that the immune status of patients with low-risk scores was better than that of patients with high-risk scores (Figure 7G). Collectively, these results suggest that the immune infiltration landscape is associated with the RM of the ICD-related genes.

Figure 7 The immune landscape between the ICD high- and low-risk groups. (A-D) Violin plots of the stromal score, immune score, ESTIMATE score, and tumor purity between the ICD high- and low-risk groups. (E) Proportion of 22 types of immune cells infiltrated between the ICD high- and low-risk groups. (F) Scores of the 22 immune cells in the high- and low-risk groups computed by the CIBERSORT algorithm. (G) Fraction of the 28 immune cells in the low- and high-risk groups computed by the ssGSEA algorithm. (H) Heatmap of the correlations among different immune cells. Blue and red illustrate positive and negative associations, respectively. ns, P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001. TME, tumor microenvironment; NK, natural killer; MDSC, myeloid-derived suppressor cell; ICD, immunogenic cell death; ssGSEA, single-sample gene-set enrichment analysis.

Finally, we also explored the relevance of the 28 different immune cells (Figure 7H). A correlation analysis was conducted to investigate the association between the prognostic ICD-related genes and the immune infiltration landscape. The results showed a remarkable association between the prognostic ICD-related genes and 22 types of immune cells as calculated by CIBERSORT; for example, FOXP3 was positively correlated with resting mast cells, and TLR4 was positively correlated with Tregs (Figure 8A). The results also revealed significant associations between the prognostic ICD-related genes. The ssGSEA, showed that 28 types of immune cells, such as LY96, were negatively correlated with macrophages (Figure 8B).

Figure 8 The association between the prognostic genes and immune cells. (A) The association between the five prognostic ICD-related genes and 22 immune cells. (B) The association between the prognostic ICD-related genes and 28 immune cells. *, P<0.05; **, P<0.01; ***, P<0.001. NK, natural killer; ICD, immunogenic cell death.

Given the significance of checkpoint inhibitor immunotherapy, we examined the expression of immune checkpoints in the two groups. Specifically, we evaluated the expression of 24 human leukocyte antigen (HLA) genes and 47 infected cell protein (ICP) genes in the two groups. The findings revealed that nearly all the ICP and HLA genes were significantly upregulated in the high-risk group, which suggested an association between this group and an immune hot phenotype. Conversely, the low-risk group was associated with an immune cold phenotype. Notably, the high-risk group had significantly elevated expression levels of molecules, such as cytotoxic T-lymphocyte-associated protein 4 (CTLA4), programmed cell death-1 (PDCD1), and Recombinant Human CD27 Ligand (CD70) (Figure 9A,9B). These findings suggest that the upregulation of immune checkpoints and immunosuppressive cytokines may suppress the immune microenvironment (IME) in high-risk patients.

Figure 9 Differential expression of HLA genes and immune checkpoints. (A,B) Box plots depicting differential expression of the immune checkpoints (A), and HLA genes (B) between the ICD high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001. HLA, human leukocyte antigen; ICD, immunogenic cell death.

Prediction of risk score and drug sensitivity

To analyze the relationship between risk score and small-molecule drug resistance due to the half-maximal inhibitory concentration (IC50) value, the pRRophetic algorithm was employed to calculate the chemotherapeutic effect of 10 commonly used small-molecule drugs (i.e., bosutinib, cisplatin, paclitaxel, cytarabine, imatinib, gefitinib, lenalidomide, pyrimethamine, rapamycin, and sunitinib) in patients with KIRC. As Figure 9 shows, the IC50s of paclitaxel, pyrimethamine, cytarabine, lenalidomide, rapamycin, cisplatin, gefitinib, sunitinib, and bosutinib were considerably more improved in the low-risk group than the high-risk group; however, the IC50 of imatinib was low in the low-risk group (Figure 10A-10J). These results illustrate the different responses to anti-tumor drugs in patients with KIRC in distinct risk subgroups and suggest the potential benefits of personalized targeted treatments for patients with KIRC.

Figure 10 Drug sensitivity analysis of the low- and high-risk groups. The IC50 values showed significant differences between the low- and high-risk groups for (A) sunitinib, (B) rapamycin, (C) pyrimethamine, (D) paclitaxel, (E) lenalidomide, (F) imatinib, (G) gefitinib, (H) cytarabine, (I) cisplatin, and (J) bosutinib. IC50, half-maximal inhibitory concentration.

As Figure 11 shows, the correlation analysis results revealed a significant association between the ICD-related genes and 12 IC50 drugs; for example, LY96 was positively correlated with cisplatin, while TLR4 was negatively correlated with paclitaxel.

Figure 11 The relationship between the expression levels of the five ICD-related genes and drug sensitivity using the ICGC database. *, P<0.05; **, P<0.01. ICD, immunogenic cell death; ICGC, International Cancer Genome Consortium.

Discussion

KIRC is the most prevalent primary malignancy in adult RCC (25). As it can be challenging to recognize the early clinical indications of KIRC, the disease is typically diagnosed in its later stages. However, patients with KIRC typically have a poor prognosis due to the limited biomarkers for early detection and prognosis prediction (26). Thus, genetic characteristics need to be discovered and predictive models developed to identify patients with distinct risks and outcomes to provide individualized treatments.

ICD has been described as a distinct type of regulated cell death that is capable of triggering a fully adaptive immune response specific to antigens via danger signals or DAMPs (27,28). The expression, function, and genetic alterations of 34 ICD-related genes were examined in the current study (19). In the tumor tissues, the majority of the differentially expressed ICD-related genes were highly expressed. Using consensus clustering, we identified two ICD subtypes based on ICD-related gene expression and found that the subtype with higher ICD expression was linked to a poor prognosis. We also identified five ICD-related genes from the 34 ICD-related prognostic genes and successfully classified the KIRC patients into high- and low-risk groups using median risk scores as the cut-off value. For the KIRC patients, this risk signature may act as an independent prognostic indicator due to its high OS predictive value. The predictive value of the model was validated by the K-M and ROC curves, which revealed that when paired with other clinical characteristics of KIRC patients, the model had the maximum prognostic predictive power. Thus, the efficiency of the nomograms in predicting the prognosis of KIRC patients was increased using a combination of risk scores and clinical characteristics. Additionally, the association between the risk score and clinicopathological features, immunological features, and drug sensitivity were investigated independently.

The five genes in the model have been implicated in tumors and other diseases. IL6 is a pro-inflammatory cytokine that has been shown to play a significant role in maintaining chronic inflammation and promoting cancer cachexia and pro-metastatic niche formation in the liver (29). FOXP3 has been reported to repress the transcription of the proto-oncogene human epidermal growth factor receptor 2 (HER2) and runt-related transcription factor 1 (RUNX1), the expression of vascular endothelial growth factor (VEGF), and the activation of the metastatic tumor antigen 1 (MTA1) anti-cancer function in breast cancer (30-35). IFNB1 has anti-proliferative and pro-apoptotic effects in Michigan Cancer Foundation-7 (MCF-7) breast cancer cells (36). LY96 interacts with TLR4 and triggers the nuclear factor kappa-B (NF-κB) pathway, which promotes the production of pro-inflammatory cytokines and adhesion molecules in colon cancer cells, thereby accelerating colon cancer growth and lung metastasis (37). The survival of colorectal cancer (CRC) patients is directly correlated with TLR4 overexpression (38-40). Experimental studies have shown that improved TLR4 activity increases CRC growth, metastasis, and immune surveillance (41,42). In clinical biopsies, higher TLR4 has been linked to the acute secretion of inflammatory cytokines, including IL6 and IL8 (43). Additionally, TLR4 facilitates the development of intestinal tumors by activating the β-catenin pathway in a manner dependent on phosphoinositide3-kinase (PI3K) (44).

The tumor IME of KIRC is characterized by substantial immune cell infiltration and various immunosuppressive effects that may severely compromise the efficacy of immunotherapy (45). We discovered 956 DEGs between the low- and high-risk groups, which helped us to better understand the gene functions and pathways in our constructed RM. According to the functional enrichment analysis, these genes were primarily associated with the immune response, which suggests that ICD is involved in TME regulation.

We examined the association between the tumor IME and risk scores, as evasion of immune destruction is a newly recognized characteristic of cancer. The high-risk patients in our study tended to have lower tumor purity and higher ESTIMATE scores, which may help to explain why these patients had a poor prognosis. The CIBERSORT algorithm findings showed that the low-risk group had higher proportions of CD4⁺ resting memory T cells, resting NK cells, monocytes, M1 macrophages, M2 macrophages, resting dendritic cells, resting mast cells, activated mast cells, and eosinophils. The high-risk group patients had higher proportions of plasma cells, CD8⁺ T cells, activated memory CD4⁺ T cells, follicular helper T cells, Tregs, and M1 macrophages. M1 macrophages have been reported to exert anti-tumor effects by producing pro-inflammatory cytokines (46). In this study, the immune infiltration analysis showed an increase in the number of infiltrating immune cells, such as CD8⁺ T, in the high-risk group. However, the high-risk score was positively correlated with programmed cell death protein 1 (PD-1) expression. Thus, while this gene can recruit immune cells into tumor tissue, the high expression of PD-1 inactivates T cells. Thus, the high-risk group continued to show an inhibitory effect in relation to the tumor immune response.

Immune checkpoint blockades (ICBs) have become a major therapeutic approach due to their ability to overturn the immunosuppressive tumor microenvironment. ICP expression is crucial for immune evasion, and ICB therapy and thus immune checkpoint inhibitors are emerging as targets following recent advances in cancer immunotherapy. In this study, most ICPs, including important ICPs, such as programmed cell death-ligand 1 (PD-L1), CD70, and CTLA4, showed substantially elevated expression levels in the highs-risk populations, indicating that high-risk patients may display higher sensitivity to ICB-based therapy. The association between ICD-related RM, immune infiltration, and ICP could be a potential research area to enhance the effectiveness of immunotherapy in solid carcinomas.

The variations in the KIRC patients’ sensitivities to the chemotherapeutic drugs may have contributed to the observed distinctions in the survival times between the high- and low-risk groups. To enhance patient prognosis and to address the fact that a lack of drug responses influences clinical decision making, we examined resistance in both risk groups and small-molecule drugs. The drug sensitivity prediction analysis results showed that the patients in the high-risk group might be more sensitive to imatinib treatment. Taken together, the results suggest that patients with high-risk scores might be more suitable candidates for immune checkpoint inhibitors and targeted therapy agents. Conversely, the selection of available treatments may be limited for patients with low-risk scores.

Conclusions

This study developed a five-gene RM and built nomograms to predict KIRC results in clinical practice. Using high- and low-risk populations, it also analyzed the IME. The correlation among immune infiltration levels, immune response levels, and pharmacologic therapy levels at both risk levels were also investigated, and a preliminary analysis of the signaling pathways involved was performed. The prognosis model constructed could better predict the survival of patients when other clinical traits, including age, gender, and stage, were taken into consideration, and the reliability of this model was verified based on information from the ICGC database. The current research offers a new and extensive viewpoint to clarify the underlying processes of KIRC prognosis and offers guidance for personalized cancer immunotherapies. Admittedly, this study had some limitations, and further research needs to be conducted to verify the results of this study.

Acknowledgments

Funding: This study was supported by the National Natural Science Foundation of China (No. 81703073).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-214/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-214/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-214/coif). Lincong Zhao is an employee of Information and Communication Branch of State Grid Hebei Electric Power Co. Ltd. However, this did not influence the manuscript. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Znaor A, Lortet-Tieulent J, Laversanne M, et al. International variations and trends in renal cell carcinoma incidence and mortality. Eur Urol 2015;67:519-30. [Crossref] [PubMed]
2. Que WC, Qiu HQ, Cheng Y, et al. PTEN in kidney cancer: A review and meta-analysis. Clin Chim Acta 2018;480:92-8. [Crossref] [PubMed]
3. Argentiero A, Solimando AG, Krebs M, et al. Anti-angiogenesis and Immunotherapy: Novel Paradigms to Envision Tailored Approaches in Renal Cell-Carcinoma. J Clin Med 2020;9:1594. [Crossref] [PubMed]
4. Linehan WM, Ricketts CJ. The Cancer Genome Atlas of renal cell carcinoma: findings and clinical implications. Nat Rev Urol 2019;16:539-52. [Crossref] [PubMed]
5. Xu WH, Xu Y, Wang J, et al. Prognostic value and immune infiltration of novel signatures in clear cell renal cell carcinoma microenvironment. Aging (Albany NY) 2019;11:6999-7020. [Crossref] [PubMed]
6. Khadirnaikar S, Kumar P, Pandi SN, et al. Immune associated LncRNAs identify novel prognostic subtypes of renal clear cell carcinoma. Mol Carcinog 2019;58:544-53. [Crossref] [PubMed]
7. Chen W, Lin W, Wu L, et al. A Novel Prognostic Predictor of Immune Microenvironment and Therapeutic Response in Kidney Renal Clear Cell Carcinoma based on Necroptosis-related Gene Signature. Int J Med Sci 2022;19:377-92. [Crossref] [PubMed]
8. Wan B, Liu B, Huang Y, et al. Prognostic value of immune-related genes in clear cell renal cell carcinoma. Aging (Albany NY) 2019;11:11474-89. [Crossref] [PubMed]
9. Galluzzi L, Vitale I, Aaronson SA, et al. Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death Differ 2018;25:486-541. [Crossref] [PubMed]
10. Zhou J, Wang G, Chen Y, et al. Immunogenic cell death in cancer therapy: Present and emerging inducers. J Cell Mol Med 2019;23:4854-65. [Crossref] [PubMed]
11. Wang X, Wu S, Liu F, et al. An Immunogenic Cell Death-Related Classification Predicts Prognosis and Response to Immunotherapy in Head and Neck Squamous Cell Carcinoma. Front Immunol 2021;12:781466. [Crossref] [PubMed]
12. Ahmed A, Tait SWG. Targeting immunogenic cell death in cancer. Mol Oncol 2020;14:2994-3006. [Crossref] [PubMed]
13. Sullivan RJ, Flaherty KT. Immunotherapy: Anti-PD-1 therapies-a new first-line option in advanced melanoma. Nat Rev Clin Oncol 2015;12:625-6. [Crossref] [PubMed]
14. Steeb T, Wessely A, Drexler K, et al. The Quality of Practice Guidelines for Melanoma: A Methodologic Appraisal with the AGREE II and AGREE-REX Instruments. Cancers (Basel) 2020;12:1613. [Crossref] [PubMed]
15. Cai J, Hu Y, Ye Z, et al. Immunogenic cell death-related risk signature predicts prognosis and characterizes the tumour microenvironment in lower-grade glioma. Front Immunol 2022;13:1011757. [Crossref] [PubMed]
16. Fu W, Ma G. Significance of immunogenic cell death-related genes in prognosis prediction and immune microenvironment landscape of patients with cutaneous melanoma. Front Genet 2022;13:988821. [Crossref] [PubMed]
17. Xu M, Lu JH, Zhong YZ, et al. Immunogenic Cell Death-Relevant Damage-Associated Molecular Patterns and Sensing Receptors in Triple-Negative Breast Cancer Molecular Subtypes and Implications for Immunotherapy. Front Oncol 2022;12:870914. [Crossref] [PubMed]
18. Kroemer G, Galluzzi L, Kepp O, et al. Immunogenic cell death in cancer therapy. Annu Rev Immunol 2013;31:51-72. [Crossref] [PubMed]
19. Garg AD, De Ruysscher D, Agostinis P. Immunological metagene signatures derived from immunogenic cancer cell death associate with improved survival of patients with lung, breast or ovarian malignancies: A large-scale meta-analysis. Oncoimmunology 2015;5:e1069938. [Crossref] [PubMed]
20. Bolstad BM, Irizarry RA, Astrand M, et al. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 2003;19:185-93. [Crossref] [PubMed]
21. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
22. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
23. Geeleher P, Cox N, Huang RS. PRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 2014;9:e107468. [Crossref] [PubMed]
24. Swets JA. Measuring the accuracy of diagnostic systems. Science 1988;240:1285-93. [Crossref] [PubMed]
25. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68:394-424. [Crossref] [PubMed]
26. Cui H, Shan H, Miao MZ, et al. Identification of the key genes and pathways involved in the tumorigenesis and prognosis of kidney renal clear cell carcinoma. Sci Rep 2020;10:4271. [Crossref] [PubMed]
27. Galluzzi L, Vitale I, Warren S, et al. Consensus guidelines for the definition, detection and interpretation of immunogenic cell death. J Immunother Cancer 2020;8:e000337. [Crossref] [PubMed]
28. Serrano-Del Valle A, Anel A, Naval J, et al. Immunogenic Cell Death and Immunotherapy of Multiple Myeloma. Front Cell Dev Biol 2019;7:50. [Crossref] [PubMed]
29. Schumacher N, Schmidt S, Schwarz J, et al. Circulating Soluble IL-6R but Not ADAM17 Activation Drives Mononuclear Cell Migration in Tissue Inflammation. J Immunol 2016;197:3705-15. [Crossref] [PubMed]
30. Lippitz BE, Harris RA. Cytokine patterns in cancer patients: A review of the correlation between interleukin 6 and prognosis. Oncoimmunology 2016;5:e1093722. [Crossref] [PubMed]
31. Lee JW, Stone ML, Porrett PM, et al. Hepatocytes direct the formation of a pro-metastatic niche in the liver. Nature 2019;567:249-52. [Crossref] [PubMed]
32. Zuo T, Wang L, Morrison C, et al. FOXP3 Is an X-Linked Breast Cancer Suppressor Gene and an Important Repressor of the HER-2/ErbB2 Oncogene. Cell 2021;184:6378. [Crossref] [PubMed]
33. Recouvreux MS, Grasso EN, Echeverria PC, et al. RUNX1 and FOXP3 interplay regulates expression of breast cancer related genes. Oncotarget 2016;7:6552-65. [Crossref] [PubMed]
34. Li X, Gao Y, Li J, et al. FOXP3 inhibits angiogenesis by downregulating VEGF in breast cancer. Cell Death Dis 2018;9:744. [Crossref] [PubMed]
35. Liu C, Han J, Li X, et al. FOXP3 Inhibits the Metastasis of Breast Cancer by Downregulating the Expression of MTA1. Front Oncol 2021;11:656190. [Crossref] [PubMed]
36. Ambjørn M, Ejlerskov P, Liu Y, et al. IFNB1/interferon-β-induced autophagy in MCF-7 breast cancer cells counteracts its proapoptotic function. Autophagy 2013;9:287-302. [Crossref] [PubMed]
37. Rajamanickam V, Yan T, Xu S, et al. Selective targeting of the TLR4 co-receptor, MD2, prevents colon cancer growth and lung metastasis. Int J Biol Sci 2020;16:1288-302. [Crossref] [PubMed]
38. Cammarota R, Bertolini V, Pennesi G, et al. The tumor microenvironment of colorectal cancer: stromal TLR-4 expression as a potential prognostic marker. J Transl Med 2010;8:112. [Crossref] [PubMed]
39. Semlali A, Reddy Parine N, Arafah M, et al. Expression and Polymorphism of Toll-Like Receptor 4 and Effect on NF-Kb Mediated Inflammation in Colon Cancer Patients. PLoS One 2016;11:e0146333. [Crossref] [PubMed]
40. Correction: Targeting immunogenic cancer cell death by photodynamic therapy: past, present and future. J Immunother Cancer 2021;9:e001926corr1.
41. Huang B, Zhao J, Li H, et al. Toll-like receptors on tumor cells facilitate evasion of immune surveillance. Cancer Res 2005;65:5009-14. [Crossref] [PubMed]
42. Fukata M, Chen A, Vamadevan AS, et al. Toll-like receptor-4 promotes the development of colitis-associated colorectal tumors. Gastroenterology 2007;133:1869-81. [Crossref] [PubMed]
43. Hu X, Fatima S, Chen M, et al. Toll-like receptor 4 is a master regulator for colorectal cancer growth under high-fat diet by programming cancer metabolism. Cell Death Dis 2021;12:791. [Crossref] [PubMed]
44. Santaolalla R, Sussman DA, Ruiz JR, et al. TLR4 activates the β-catenin pathway to cause intestinal neoplasia. PLoS One 2013;8:e63298. [Crossref] [PubMed]
45. Díaz-Montero CM, Rini BI, Finke JH. The immunology of renal cell carcinoma. Nat Rev Nephrol 2020;16:721-35. [Crossref] [PubMed]
46. Wculek SK, Dunphy G, Heras-Murillo I, et al. Metabolism of tissue macrophages in homeostasis and pathology. Cell Mol Immunol 2022;19:384-408. [Crossref] [PubMed]