Research Paper Volume 15, Issue 12 pp 5304—5338

Construction and validation of a model based on immunogenic cell death-associated lncRNAs to predict prognosis and direct therapy for kidney renal clear cell carcinoma

Chenxi Cai1,2, , Kexin Shu2,3, , Wanying Chen2,3, , Jiatong Ding2,3, , Zishun Guo1, , Yiping Wei1, &, , Wenxiong Zhang1, &, ,

  • 1 Department of Thoracic Surgery, The Second Affiliated Hospital of Nanchang University, Nanchang 330006, China
  • 2 Jiangxi Medical College, Nanchang University, Nanchang 330006, China
  • 3 Department of Urinary Surgery, The Second Affiliated Hospital of Nanchang University, Nanchang 330006, China

Received: January 31, 2023       Accepted: May 9, 2023       Published: June 27, 2023
How to Cite

Copyright: © 2023 Cai et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Background: Immunogenic cell death (ICD) is an important part of the antitumor effect, yet the role played by long noncoding RNAs (lncRNAs) remains unclear. We explored the value of ICD-related lncRNAs in tumor prognosis assessment in kidney renal clear cell carcinoma (KIRC) patients to provide a basis for answering the above questions.

Methods: Data on KIRC patients were obtained from The Cancer Genome Atlas (TCGA) database, prognostic markers were identified, and their accuracy was verified. An application-validated nomogram was developed based on this information. Furthermore, we performed enrichment analysis, tumor mutational burden (TMB) analysis, tumor microenvironment (TME) analysis, and drug sensitivity prediction to explore the mechanism of action and clinical application value of the model. RT-qPCR was performed to detect the expression of lncRNAs.

Results: The risk assessment model constructed using eight ICD-related lncRNAs provided insight into patient prognoses. Kaplan-Meier (K-M) survival curves showed a more unfavorable outcome in high-risk patients (p<0.001). The model had good predictive value for different clinical subgroups, and the nomogram constructed based on this model worked well (risk score AUC=0.765). Enrichment analysis revealed that mitochondrial function-related pathways were enriched in the low-risk group. The adverse prognosis of the higher-risk cohort might correspond to a higher TMB. The TME analysis revealed a higher resistance to immunotherapy in the increased-risk subgroup. Drug sensitivity analysis can guide the selection and application of antitumor drugs in different risk groups.

Conclusions: This prognostic signature based on eight ICD-associated lncRNAs has significant implications for prognostic assessment and treatment selection in KIRC.


Kidney renal clear cell carcinoma (KIRC) is the most common subtype of renal cell carcinoma, constituting approximately 75%-85% of renal cell carcinoma cases [1]. Its morbidity is rising each year, and patients often present with postoperative metastases, which are associated with great patient hardship [2]. The American Joint Committee on Cancer (AJCC) tumor node metastasis (TNM) classification system is often used for staging and is used to divide patients into stages I, II, III, and IV for prognostic assessment [3]. However, the predictions of the prognosis of patients treated based on this system are not accurate [4]. Therefore, novel prognostic models are urgently needed to guide the treatment of patients with KIRC. Biomarker-based prognostic models have shown great potential in recent years for tumor patient prognostic assessment.

Immunogenic cell death (ICD), one of the important modalities of regulatory cell death, promotes antitumor effects by triggering the activation of cytotoxic T cells [5]. Simultaneously, long noncoding RNAs (lncRNAs) are indispensable for performing this function, and lncRNA-related models have been successfully constructed for numerous cancer types. lncRNA models constructed by Liu Z et al. were validated in the prognostic assessment of colorectal cancer patients [6]. lncRNA contributions in tumor therapy were reviewed by Eptaminitaki GC et al. [7]. Liang YL et al. constructed a tumor immune heterogeneity-associated lncRNA prognostic model to determine the long-term prognosis of patients with nasopharyngeal carcinoma [8]. However, a KIRC prognostic assessment model based on ICD-related lncRNAs is still not available.

Thus, we constructed an ICD-associated lncRNA-based prognostic model for assessing the prognosis of patients with KIRC and appraised its clinical application value through enrichment analysis, tumor mutational burden (TMB) analysis, tumor microenvironment (TME) differential analysis, and drug sensitivity prediction.

Materials and Methods

Data sources and access

The University of California, Santa Cruz (UCSC) Xena database (, until November 1st, 2022) was the primary source of data for this study [9]. The Cancer Genome Atlas (TCGA) database in Xena provided patient information, including transcriptomic data for 607 KIRC cases and clinicopathological information for 979 KIRC clinical patients. Tumor somatic cell mutation data were obtained from the TCGA database (, until November 1st, 2022). Transcriptomic data were obtained in two data formats, HTSeq-Counts and HTSeq-FPKM, which are used for different types of data analysis. We selected patients with complete transcriptional and clinical information and ensured that they were identified in both datasets. Ultimately, a total of 597 samples were analyzed in this study.

Select ICD-associated genes and ICD-related lncRNAs

ICD-associated genes were identified from the GeneCards website (, until November 1st, 2022). We conducted a search using immunogenic cell death as a keyword, and target genes were selected for subsequent analysis [10]. In addition, we conducted Pearson correlation analysis to investigate the associations among ICD-associated genes and all lncRNAs and acquired ICD-related lncRNAs (filter conditions indicated as |correlation coefficient| > 0.3 and p < 0.05). We ran variance analysis on the results obtained from Pearson correlation analysis using the “DESeq2” package [11], with the retention of lncRNAs with a p value < 0.05 and |log2-fold change > 3|.

All KIRC patients were randomly assigned 1:1 to the training cohort, which was used for model construction, or the test cohort, which was used to verify the model. Then, ICD-associated lncRNAs relevant to patient survival were short-listed in the training cohort using univariate Cox regression analysis. We used least absolute shrinkage and selection operator (LASSO) regression analysis to avoid overfitting. Finally, ICD-related lncRNAs were identified through multivariate Cox regression analysis.

KIRC ICD-associated lncRNA prognostic model

We used the results from predictive model building. The format of the risk score was as follows: risk score = coefficient (lncRNA1) × expression (lncRNA1) + coefficient (lncRNA2) × expression (lncRNA2) + coefficient (lncRNA3) × expression (lncRNA3) +......+ coefficient (lncRNAn) × expression (lncRNAn). Using the midpoint of the risk score as its threshold, the training group, the test group and all patients were sorted into high- and low-risk cohorts for survival analysis. Kaplan-Meier (K-M) survival analysis utilizing “survival” and “survminer” was conducted to examine the variation in overall survival (OS) between risk categories in the training cohort, test cohort, and total cohort. Receiver operating characteristic (ROC) values at one, three, and five years were extrapolated to evaluate the predictive performance of the signature.

Principal component analysis (PCA) using the expression profiles of all genes, all lncRNAs, lncRNAs associated with ICD, and lncRNAs in the selected prognostic models was conducted to validate the subgrouping effect. In addition, univariate and multifactor Cox independent prognostic analyses were conducted to assess risk scores and clinical data to validate the predictive value of the risk model. After that, K-M survival analysis was conducted to probe the model feasibility under dissimilar clinical characteristics.

Creation of a nomogram for predicting patient survival

We used the R programming language “RMS” to create a prognostic nomogram using age, risk group, and tumor stage as factors to predict patient outcomes at 1, 3, and 5 years [12]. We also determined the value of this nomogram in prognosis prediction.

Enrichment analysis

KEGG enrichment analysis was applied to optimize the functions that might be interrelated with the genetic set. We divided the cases into high- and low-risk groups under the simulations constructed from 8 ICD-related lncRNAs, and enrichment analysis was used to filter the pathways with a p value<0.05 and false discovery rate (FDR)<0.25 using GSEA_4.3.2 software (downloaded at [13].

Tumor mutational burden contributes to the prognostic evaluation of tumors

We obtained tumor somatic mutation data from the TCGA database. We manipulated the data using the “TCGAbiolinks” package, sketched waterfall plots using the R programming language “maftools”, and calculated TMB [14, 15].

Evaluation of TME and immuno-infiltration

The R software package “ESTIMATE” was used as an analysis vehicle to calculate any difference in stromal score, ESTIMATE score, immune score, and tumor purity between the high- and low-risk subgroups. We obtained further information through the Tumor Immune Estimation Resource (TIMER) 2.0 website (, until November 1st, 2022) using “TIMER”, “CIBERSORT”, “CIBERSORT-ABS”, “QUANTISEQ”, “XCELL”, and “EPIC”, together with “MCPCOUNTER”, which are six methods used to optimize the correlation between individual immune cell types and risk scores [16]. To detail the constitutive shifts in immune cells, we undertook an immune infiltration analysis of 22 immune cell lineages.

ICD-related signature in immunotherapy versus chemotherapy

To fully address the distinctions between KIRC patients in different risk groups presenting with tumor immune dysfunction and rejection, the Tumor Immune Dysfunction and Exclusion (TIDE) website (, until November 1st, 2022) was used to retrieve the TIDE scores and related data of KIRC patients [17]. Therefore, we conducted single-sample gene set enrichment analysis (ssGSEA) using the “GSVA” R package. Half-maximal inhibitory concentrations (IC50) of typical antineoplastic drugs were determined using the R software “oncoPredict” package [18]. The IC50 values were compared between different classes utilizing a Wilcoxon signed-rank test.

Validation by RT-qPCR

Twenty-one KIRC samples were obtained from the Second Affiliated Hospital of Nanchang University, and patients were divided into high- and low-risk groups. The study was approved by the Ethics Committee of the Second Affiliated Hospital of Nanchang University, and all participants gave their informed consent. We extracted RNA from each KIRC tissue sample using TRIzol reagent (Life Technologies CA, USA) and randomly selected samples for RT-qPCR analysis. The experiments were performed using BlazeTaq SYBR Green qPCR master mix (GeneCopoeia, Guangzhou, China) and the Applied Biosystems 7500 Fast Real-Time PCR System (Applied Biosystems). All RNAs of every sample were analyzed in three independent experiments. The primers for ICD-associated lncRNAs are shown in Supplementary Table 1. The relative expression of lncRNAs was calculated using the 2^(-ΔΔCt) method.

We used the Human Protein Atlas ( database to compare the protein expression levels of selected ICD-related genes in KIRC tissues with those in normal tissues.

Availability of data and material

The data sets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.


Screening for ICD-associated lncRNAs in KIRC patients

The overall study design of this paper is displayed in a flow chart (Figure 1). First, we identified ICD-related genes and constructed a protein–protein interaction (PPI) network associated with 16 of these genes by using the STRING website (Figure 2A). Further Pearson correlation analysis was executed, from which we derived 9096 ICD-associated lncRNAs. Next, we narrowed the pool down to 547 lncRNAs that were shown by differential analysis to be differentially expressed in tumor tissues (Figure 2B).

Flow chart.

Figure 1. Flow chart.

Filter for ICD-related lncRNAs. PPI network among 16 ICD-related genes (A); Volcano plot of differentially expressed ICD-associated lncRNAs (B); LASSO regression analysis (C, D).

Figure 2. Filter for ICD-related lncRNAs. PPI network among 16 ICD-related genes (A); Volcano plot of differentially expressed ICD-associated lncRNAs (B); LASSO regression analysis (C, D).

Based on difference analysis, we rerandomized 526 patients into training and test groups (Table 1). In the training group, a single-factor Cox model was used to screen 180 lncRNAs associated with differential prognosis using p<0.05 as the cutoff value (Supplementary Table 2). Follow-up LASSO analysis revealed 14 differentially expressed ICD-related lncRNAs (Figure 2C, 2D) (Supplementary Table 3), and multifactorial Cox regression analysis was used to further screen the 8 most significant lncRNAs (Supplementary Figure 1A). We subsequently explored the correlations among these 8 lncRNAs (Supplementary Figure 1B) and their correlation with ICD-related genes (Supplementary Figure 1C) that were found to exhibit a marked difference in up- and downregulation in tumor tissues (Supplementary Figure 1D).

Table 1. Clinical information of the patients in the test and training groups.

CharacteristicsTrain cohort (n=264)Test cohort (n=262)Entire cohort (n=526)
Stage I12848.4813350.7626149.62
Stage II3011.362710.315710.84
Stage III6223.486123.2812323.38
Stage IV4215.914015.278215.59
T stage
M stage
N stage
Black or African American259.472911.075410.27
Abbreviation: T stage, Tumor stage; N stage, Node stage; M stage, metastasis stage.

Building and validating the prognostic model

Based on the 8 ICD-related lncRNAs, we calculated risk scores for different patients: risk score= coefficient(AP000439.3)×expression (AP000439.3)+ coefficient (RP11.1151B14.5)×expression(RP11.1151B14.5)+coefficient(RP11.479J7.2)×expression(RP11.479J7.2)+coefficient(AC099552.4)×expression(AC099552.4)+coefficient(RP11.19E11.1)×expression(RP11.19E11.1)+coefficient(CTB.33O18.1)×expression(CTB.33O18.1)+coefficient(RP11.339D23.1)×expression(RP11.339D23.1)+coefficient(LINC01192)×expression(LINC01192). Clinical information for all patients in the high- and low-risk clusters is presented in Table 2. In the training cohort, the test cohort, and the total cohort, survival was much worse in the high-risk group than in the low-risk group (all p values<0.001) (Figure 3A3C), and the time-dependent ROC curve demonstrated that the model’s predictive value was excellent (AUCs were 0.831, 0.775, and 0.796 for the training cohort at one, three, and five years, respectively; 0.710, 0.688, and 0.768 for the test cohort at one, three, and five years, respectively; and 0.775, 0.743, and 0.786 for the total cohort at one, three, and five years, respectively) (Figure 3D3F). In different cohorts, the manifestations of 8 lncRNAs in patients with different risks (Supplementary Figure 2A2C), the risk curves (Supplementary Figure 2D2F), the risk distribution plot (Supplementary Figure 2G2I), and the scatter plot (Supplementary Figure 2J2L) showed the forecasting value of the model.

Table 2. Clinical information for 526 patients in different risk categories.

CharacteristicsHigh-risk group (n=263)Low-risk group (n=263)
Stage I10038.0216161.22
Stage II2710.273011.41
Stage III7428.144918.62
Stage IV6022.81228.37
T stage
M stage
N stage
Black or African American2911.03259.51
Abbreviation: T stage, Tumor stage; N stage, Node stage; M stage, metastasis stage.
The model prediction effect is validated by the training group, test group, and entire group. K-M analysis (A–C) and Time-dependent ROC curves (D–F) to compare the survival of the high-risk group and low-risk group.

Figure 3. The model prediction effect is validated by the training group, test group, and entire group. K-M analysis (AC) and Time-dependent ROC curves (DF) to compare the survival of the high-risk group and low-risk group.

The outcomes of principal component analysis (PCA) (Figure 4A4D) indicated that the high- and low-risk subgroups showed obvious categorical clustering. Univariate Cox independent prognostic analysis indicated that age (HR=1.68, 95% CL=1.24-2.28, p<0.001), tumor stage (HR=1.92, 95% CL=1.68-2.20, p<0.001), and risk score (HR=1.44, 95% CL=1.34-1.54, p<0.001) were standalone risk factors (Figure 4E). Likewise, multifactorial independent prognostic analysis showed that age (HR=1.57, 95% CL=1.15-2.15, p=0.004), tumor stage (HR=1.73, 95% CL=1.50-1.99, p<0.001), and risk score (HR=1.34, 95% CL=1.25-1.44, p<0.001) were standalone risk factors (Figure 4F).

PCA and independent prognostic analysis of the signature. PCA based on all genes (A), all lncRNAs (B), ICD-related lncRNAs (C), and risk signature (D); Univariate (E) and multivariate (F) independent prognostic analysis.

Figure 4. PCA and independent prognostic analysis of the signature. PCA based on all genes (A), all lncRNAs (B), ICD-related lncRNAs (C), and risk signature (D); Univariate (E) and multivariate (F) independent prognostic analysis.

The expression heatmap containing different clinical data, risk groupings, and prognostic model-related lncRNAs (Supplementary Figure 3A) and the survival curves of patients in a variety of clinical states (Figure 5) illustrate the significance of the model predictions (stage I-II: p<0.001, stage III-IV: p<0.001; female: p<0.001, male: p<0.001; age<=65: p<0.001, age>65: p=0.002).

Further validation of model effects. Survival curves of patients in different clinical states (A–F).

Figure 5. Further validation of model effects. Survival curves of patients in different clinical states (AF).

Clinical OS prediction nomogram

The scatter plot showed a positive correlation with respect to tumor stage and risk scores (p<0.001) (Supplementary Figure 3C). Age >65 years also tended to be a factor for an increased risk score (p=0.096) (Supplementary Figure 3B). Decision curve analysis (DCA) (Figure 6A) and ROC curves (Figure 6B) (risk AUC=0.765, age AUC=0.646, sex AUC=0.500, and stage AUC=0.815) showed that risk had a superior prediction value compared to most clinical information. A prospective estimator of KIRC patients based on age, risk, sex, and disease staging was conceived (Figure 6C), patients were used to verify its effectiveness, and the results indicated good performance (Supplementary Figure 4). The forecasting value of the nomogram combined with the risk model (Supplementary Figure 4D) was higher than that of the nomogram without the risk model (Supplementary Figure 4E).

Nomogram predicts patient prognosis. Decision curve to test for forecast value (A); ROC curves containing different clinical information (B); A clinical prognosis nomogram is constructed by age, gender, risk, and stage together (C). Nomogram with (D) and without (E) risk model.

Figure 6. Nomogram predicts patient prognosis. Decision curve to test for forecast value (A); ROC curves containing different clinical information (B); A clinical prognosis nomogram is constructed by age, gender, risk, and stage together (C). Nomogram with (D) and without (E) risk model.

Enrichment analysis

KEGG enrichment analysis revealed the related functions of variably expressed genes, including ubiquitin-mediated proteolysis, glycolipid biosynthesis, galactose metabolism, ethoxylate, and dicarboxylate metabolism (Supplementary Figure 5). GSEA demonstrated that different pathways were enriched in different gene sets. Pathway functions enriched in the low-risk category in the GOBP library included mitochondrial gene expression, assembly of the mitochondrial respiratory chain complex, mitochondrial translation, neurotransmitter reuptake, and mitochondrial electron transport from NADH to ubiquinone (Supplementary Figure 6). All the pathway information is shown in Supplementary Table 4.

Tumor mutational burden

The TMB in patients was determined and found to be higher in the high-risk subgroup (Figure 7A). Based on this, we generated a waterfall plot of the top 20 mutant genes according to different subclusters (Figure 7B, 7C). The five most commonly mutated genes in the high-risk patients were VHL (45%), PBRM1 (38%), TTN (17%), BAP1 (16%), and SETD2 (16%). VHL (48%), PBRM1 (42%), TTN (15%), SETD2 (8%) and MUC16 (7%) were prone to mutation in low-risk patients.

Tumor mutation burden in different risk groups. Percentage bar graph showing TMB for different risk subgroups (A); High-risk group waterfall chart (B); Low-risk group waterfall chart (C).

Figure 7. Tumor mutation burden in different risk groups. Percentage bar graph showing TMB for different risk subgroups (A); High-risk group waterfall chart (B); Low-risk group waterfall chart (C).

Tumor immune infiltration status analysis

The TME analysis suggested that the immune score (Figure 8A) and the ESTIMATE score (Figure 8C) were higher in the high-risk segment (p<0.01), but the stromal score (Figure 8B) did not show dramatic differences. Patients at high risk had significantly lower tumor purity (p<0.01) (Figure 8D). We hypothesized that an immunosuppressive microenvironment was present in the high-risk subgroup that weakened antitumor immunity. Additionally, we showed an association between immune cells and the risk score under different algorithms (Figure 8E). Naive B cells (p<0.05), plasma cells (p<0.05), follicular helper T cells (p<0.01), regulatory T cells (Tregs) (p<0.0001), and M0 macrophages (p<0.01) constituted a greater proportion within the high-risk group, and resting memory CD4 T cells (p<0.05), monocytes (p<0.01), M1 macrophages (p<0.0001), and resting mast cells (p< 0.0001) showed comparatively higher expression in the low-risk group (Figure 8F). The relationships among immune cells and risk scores are shown in Supplementary Figure 7. Accordingly, APC coinhibition (p<0.01) and type II IFN response (p<0.001) functions were inhibited in high-risk KIRC patients, while parainflammation (p<0.01) and T-cell costimulation functions were improved (Figure 8G). It can be speculated that tumor development in high-risk KIRC patients is facilitated by both T-cell parainflammation and costimulation.

Analysis of tumor immune microenvironment. Violin plots of differences in immune scores (A), stromal scores (B), ESTIMATE scores (C), and tumor purity (D) for different risk subgroups; Bubble plots of correlations between immune cells and risk scores under six algorithms (E); Proportions of 22 immune cells in two subgroups under the CIBERSORT algorithm (F); single sample gene set enrichment analysis (G). *p

Figure 8. Analysis of tumor immune microenvironment. Violin plots of differences in immune scores (A), stromal scores (B), ESTIMATE scores (C), and tumor purity (D) for different risk subgroups; Bubble plots of correlations between immune cells and risk scores under six algorithms (E); Proportions of 22 immune cells in two subgroups under the CIBERSORT algorithm (F); single sample gene set enrichment analysis (G). *p < 0.05, **p < 0.01, ***p < 0.001.

Benefits of promotional models in the treatment of KIRC

T-cell exclusion and T-cell dysfunction, together with TIDE scores, were consistently higher in the high-risk patients with ICD-related signatures (Supplementary Figure 8). This may indicate that there is more potential for tumors in the high-risk subgroup to exhibit immune escape, complicating treatment. The antitumor drugs with different mechanisms of action to which tumors in the low-risk subgroup are sensitive are shown in Supplementary Table 5, those to which tumors in the higher-risk subgroup are sensitive are shown in Supplementary Table 6, and those without significant intergroup susceptibility differences are shown in Supplementary Table 7. The high-risk group appears to be more sensitive to drugs that act on the PI3K/mTOR signaling pathway and metabolism pathway, which can guide drug selection.

In vitro experimental validation of risk models

The protein expression levels of selected ICD-related genes in KIRC tissues and normal tissues were visualized with immunohistochemical staining images in the HPA database (Supplementary Figure 9A).

RT-qPCR results showed that among the eight ICD-related lncRNAs, AP000439.3 was highly expressed in the low-risk group, and the remaining seven lncRNAs were highly expressed in the high-risk group (Supplementary Figure 9B).


KIRC is one of the major pathological types of renal cell carcinoma, which often spreads and metastasizes [1]. Nevertheless, the etiology and pathogenesis of this tumor remain to be further explored [19]. Hence, it is particularly important to accurately assess a patient’s prognosis and provide appropriate treatment. However, clinical TNM staging is not able to accurately assess the prognosis of patients with KIRC [4]. Accordingly, there is a great need to construct a new prognostic assessment model. ICD is an important component of regulated cell death and is involved in the antitumor process [5]. Thus, we screened ICD-related lncRNAs for the construction of a KIRC prognostic assessment model and explored their possible molecular mechanisms and clinical applications. The outcome indicated that the risk score derived from 8 ICD-associated lncRNAs could be used as a standalone predictive factor and that patients in the high-risk subgroup had a worse prognosis. The nomogram constructed in accordance with this had good predictive value. Enrichment analysis showed that mitochondria-associated pathways might be relevant in the low-risk subgroup. The model provides a reference for antitumor drug selection for KIRC patients.

ICD involves the release of danger-associated molecular patterns (DAMPs) from apoptotic tumor cells, which activate immune cells, thus promoting the antitumor effect of immune cells [5]. ICD-based prognostic models have good predictive value in other cancer types. Cai J used an ICD-based evaluation model to verify the prognosis of low-grade glioma patients with good results [20], and an ICD-related prognostic model constructed by Ma J. and team to forecast the prognosis of patients suffering from hepatocellular carcinoma also achieved satisfactory results [21]. We downloaded sample data from TCGA-KIRC patients, performed Pearson correlation analysis to identify ICD-related genes, and then carried out difference analysis to select test groups from the obtained tumor samples to conduct subsequent one-way Cox regression analysis, LASSO regression analysis, and multifactor regression analysis to obtain ICD-related lncRNAs to construct prognostic models. A nomogram was constructed to forecast the 1-, 3-, and 5-year survival rates of patients. The survival rates of patients in the ICD-related high-risk subgroup according to the prognostic model were all worse than those of patients in the ICD-related low-risk subgroup, while the risk scores had a better predictive effect than traditional tumor staging, which we speculated was due to the worse ICD effect in the upper-risk subgroup, leading to a worse prognosis. When the risk scores were applied, the nomogram predictions were better. Among ICD-related lncRNAs, AP000439.3 is regulated by estrogen receptor (ER) and can regulate CCND1 expression through enhancement of estrogen receptors, thereby inhibiting cell cycle progression and cell proliferation [22]. In contrast, LINC01192 expression is upregulated in triple-negative breast cancer and is associated with the low survival likelihood of patients with triple-negative breast cancer [23]. This may give better predictive value to the risk score. Overall, our ICD-associated lncRNA prognostic model showed good predictive properties in KIRC patients and had more potential than conventional assessment methods.

We carried out GSEA and TMB analysis to detect the possible mechanisms involved. GSEA revealed numerous pathways related to mitochondrial function enriched in the low-risk subgroup. It has been shown that the enzyme RIPK3 can mediate signaling between mitochondria and the immune system to initiate antitumor immunity [24]. Impaired oxidative phosphorylation due to mitochondrial defects leads to cellular carcinogenesis [25]. More vigorous aerobic glycolysis in cancerous tissues provides energy to the tissue [26]. Additionally, a higher TMB usually results in a worse prognosis for patients with KIRC [27].

We further explored the contribution of this prognostic model to clinical drug selection. It has been shown that plasma cells can produce immunoglobulins and inhibit cell growth in the early stages of disease. At an early stage, pathological IgG can enter tumor cells through the AP2 complex and degrade overexpressed proteins through the TRIM21-mediated ubiquitin pathway, thus achieving antitumor effects [28]. M1 isoforms of tumor-associated macrophages (TAMs) can enhance antitumor immunity, and mast cell resting tends to lead to a better prognosis in KIRC patients [29, 30]. The higher content of the above cells in the low-risk subgroup may prolong patient survival. Tregs can use CTLA4 to inhibit the costimulatory signaling molecules CD80 and CD86, secrete suppressive cytokines, and directly kill effector T cells, creating an immunosuppressive microenvironment [31]. High Treg expression in high-risk groups may cause worse patient prognosis. All these immune cells point the way to immunotherapy. Additionally, the high-risk subgroup had higher sensitivity to drugs targeting the PI3K/mTOR pathway. Cellular responses dominated by the PI3K/AKT/mTOR pathway are frequently seen in KIRC and are associated with tumor progression [32]. This might be connected to the higher sensitivity of the high-risk group. Therefore, this model is expected to be useful in the selection of antineoplastic drugs for KIRC patients.

We consider the following advantages of our study. We did not find other studies reporting an ICD-associated lncRNA model that had been successfully developed to predict the outcome of KIRC patients. Nevertheless, to verify the validity of our model, we compared it with other similar studies. First, the autophagy-related lncRNA model constructed by Xuan Y et al. showed predictive value in KIRC patients, but we further calculated the tumor mutational load for our constructed model and explored its use in applications such as antitumor drug selection [33]. Second, compared with the focal death-related lncRNA model constructed by Zhou X et al., we validated our model more comprehensively using all patients, and the results were more convincing [34]. Third, the prognostic model we constructed (AUC=0.765) was more meaningful than the metastasis-related lncRNA prognostic model constructed by Dou Q et al. (AUC=0.755) for the prognostic assessment of KIRC [35]. We must acknowledge that our study still has shortcomings. The information in this study was obtained from one of many databases and was not validated with external data. Thus, further clinical experiments are needed to validate this study in the future.


Overall, we developed a prognostic model with eight ICD-associated lncRNAs and constructed a nomogram, which was shown to be a valuable guide in prognostic assessment of and treatment selection for KIRC. The inferior prognosis in the high-risk cohort may be correlated with mitochondria-associated pathways and higher TMB. Due to various shortcomings, this study awaits further basic research to explore the relevant mechanisms and clinical control studies to clarify the value of the model in drug selection. Furthermore, the results of the prognostic model based on clinical samples need to be validated in further trials.


ATP: Adenosine triphosphate; CCND1: Recombinant Cyclin D1; DCA: Decision curve analysis; FDR: False discovery rate; GSEA: Gene set enrichment analysis; KEGG: Kyoto Encyclopedia of Genes and Genomes; H: high; HR: Hazard ratios; IC50: half maximal inhibitory concentration; ICD: immunogenic cell death; KIRC: kidney renal clear cell carcinoma; K-M survival analysis: Kaplan-Meier survival analysis; L: low; LASSO: Least absolute shrinkage and selection operator; lncRNAs: Long non-coding RNA; OS: Overall survival; PCA: Principal component analysis; PPI: protein-protein interaction; ROC: receiver operating characteristic; TCGA: The Cancer Genome Atlas; TMB: Tumor Mutation Burden; TME: tumor microenvironment; TAM: tumor-associated macrophage.

Author Contributions

Wenxiong Zhang had full access to all the data in the manuscript and takes responsibility for the integrity of the data and the accuracy of the data analysis. All authors read and approved the final manuscript.


The authors thank professor Jianjun Xu, MD (Department of Thoracic Surgery, The Second Affiliated Hospital of Nanchang University) for his statistical advice. The authors also thank the support of the National Natural Science Foundation of China and Natural Science Foundation of Jiangxi Province.

Conflicts of Interest

The authors certify that there are no conflicts of interest regarding this manuscript.

Ethical Statement and Consent

The current study investigated the publicly available data. The study was approved by the Ethics Committee of the Second Affiliated Hospital of Nanchang University, and all participants gave their informed consent to use KIRC samples. All methods were carried out in accordance with the Declaration of Helsinki.


This study was supported by National Natural Science Foundation of China (NSFC, Grant number: 81560345) and Natural Science Foundation of Jiangxi Province (Grant number: 20212BAB206050). Role of the Funding: The funding had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Editorial Note


This corresponding author has a verified history of publications using a personal email address for correspondence.


  • 1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin. 2022; 72:7–33. [PubMed]
  • 2. Vuong L, Kotecha RR, Voss MH, Hakimi AA. Tumor Microenvironment Dynamics in Clear-Cell Renal Cell Carcinoma. Cancer Discov. 2019; 9:1349–57. [PubMed]
  • 3. Motzer RJ, Jonasch E, Agarwal N, Alva A, Baine M, Beckermann K, Carlo MI, Choueiri TK, Costello BA, Derweesh IH, Desai A, Ged Y, George S, et al. Kidney Cancer, Version 3.2022, NCCN Clinical Practice Guidelines in Oncology. J Natl Compr Canc Netw. 2022; 20:71–90. [PubMed]
  • 4. Klatte T, Rossi SH, Stewart GD. Prognostic factors and prognostic models for renal cell carcinoma: a literature review. World J Urol. 2018; 36:1943–52. [PubMed]
  • 5. Li X, Liu Z, Zhang A, Han C, Shen A, Jiang L, Boothman DA, Qiao J, Wang Y, Huang X, Fu YX. NQO1 targeting prodrug triggers innate sensing to overcome checkpoint blockade resistance. Nat Commun. 2019; 10:3251. [PubMed]
  • 6. Liu Z, Liu L, Weng S, Guo C, Dang Q, Xu H, Wang L, Lu T, Zhang Y, Sun Z, Han X. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. Nat Commun. 2022; 13:816. [PubMed]
  • 7. Eptaminitaki GC, Stellas D, Bonavida B, Baritaki S. Long non-coding RNAs (lncRNAs) signaling in cancer chemoresistance: From prediction to druggability. Drug Resist Updat. 2022; 65:100866. [PubMed]
  • 8. Liang YL, Zhang Y, Tan XR, Qiao H, Liu SR, Tang LL, Mao YP, Chen L, Li WF, Zhou GQ, Zhao Y, Li JY, Li Q, et al. A lncRNA signature associated with tumor immune heterogeneity predicts distant metastasis in locoregionally advanced nasopharyngeal carcinoma. Nat Commun. 2022; 13:2996. [PubMed]
  • 9. Goldman MJ, Craft B, Hastie M, Repečka K, McDade F, Kamath A, Banerjee A, Luo Y, Rogers D, Brooks AN, Zhu J, Haussler D. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020; 38:675–8. [PubMed]
  • 10. Xu C, Li Y, Su W, Wang Z, Ma Z, Zhou L, Zhou Y, Chen J, Jiang M, Liu M. Identification of immune subtypes to guide immunotherapy and targeted therapy in clear cell renal cell carcinoma. Aging (Albany NY). 2022; 14:6917–35. [PubMed]
  • 11. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15:550. [PubMed]
  • 12. Chen Q, Wang S, Lang JH. Development and validation of nomogram with tumor microenvironment-related genes and clinical factors for predicting overall survival of endometrial cancer. J Cancer. 2021; 12:3530–8. [PubMed]
  • 13. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. [PubMed]
  • 14. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018; 28:1747–56. [PubMed]
  • 15. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, Ceccarelli M, Bontempi G, Noushmehr H. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016; 44:e71. [PubMed]
  • 16. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020; 48:W509–14. [PubMed]
  • 17. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24:1550–8. [PubMed]
  • 18. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021; 22:bbab260. [PubMed]
  • 19. Huang Q, Zi H, Luo L, Li X, Zhu C, Zeng X. Secular trends of morbidity and mortality of prostate, bladder, and kidney cancers in China, 1990 to 2019 and their predictions to 2030. BMC Cancer. 2022; 22:1164. [PubMed]
  • 20. Cai J, Hu Y, Ye Z, Ye L, Gao L, Wang Y, Sun Q, Tong S, Yang J, Chen Q. Immunogenic cell death-related risk signature predicts prognosis and characterizes the tumour microenvironment in lower-grade glioma. Front Immunol. 2022; 13:1011757. [PubMed]
  • 21. Ma J, Kuang L, Zhao R. Establishing a signature based on immunogenic cell death-related gene pairs to predict immunotherapy and survival outcomes of patients with hepatocellular carcinoma. Aging (Albany NY). 2022; 14:9699–714. [PubMed]
  • 22. Zhang Y, Wang DL, Yan HY, Liao JY, He JH, Hu KS, Deng WX, Wang YJ, Xing HT, Koeffler HP, Yin D. Genome-wide study of ER-regulated lncRNAs shows AP000439.3 may function as a key regulator of cell cycle in breast cancer. Oncol Rep. 2017; 38:3227–37. [PubMed]
  • 23. Li XX, Wang LJ, Hou J, Liu HY, Wang R, Wang C, Xie WH. Identification of Long Noncoding RNAs as Predictors of Survival in Triple-Negative Breast Cancer Based on Network Analysis. Biomed Res Int. 2020; 2020:8970340. [PubMed]
  • 24. Kang YJ, Bang BR, Han KH, Hong L, Shim EJ, Ma J, Lerner RA, Otsuka M. Regulation of NKT cell-mediated immune responses to tumours and liver inflammation by mitochondrial PGAM5-Drp1 signalling. Nat Commun. 2015; 6:8371. [PubMed]
  • 25. Srinivasan S, Guha M, Dong DW, Whelan KA, Ruthel G, Uchikado Y, Natsugoe S, Nakagawa H, Avadhani NG. Disruption of cytochrome c oxidase function induces the Warburg effect and metabolic reprogramming. Oncogene. 2016; 35:1585–95. [PubMed]
  • 26. Hensley CT, Faubert B, Yuan Q, Lev-Cohain N, Jin E, Kim J, Jiang L, Ko B, Skelton R, Loudat L, Wodzak M, Klimko C, McMillan E, et al. Metabolic Heterogeneity in Human Lung Tumors. Cell. 2016; 164:681–94. [PubMed]
  • 27. Valero C, Lee M, Hoen D, Wang J, Nadeem Z, Patel N, Postow MA, Shoushtari AN, Plitas G, Balachandran VP, Smith JJ, Crago AM, Long Roche KC, et al. The association between tumor mutational burden and prognosis is dependent on treatment context. Nat Genet. 2021; 53:11–5. [PubMed]
  • 28. Chen J, Tan Y, Sun F, Hou L, Zhang C, Ge T, Yu H, Wu C, Zhu Y, Duan L, Wu L, Song N, Zhang L, et al. Single-cell transcriptome and antigen-immunoglobin analysis reveals the diversity of B cells in non-small cell lung cancer. Genome Biol. 2020; 21:152. [PubMed]
  • 29. Cordova AF, Ritchie C, Böhnert V, Li L. Human SLC46A2 Is the Dominant cGAMP Importer in Extracellular cGAMP-Sensing Macrophages and Monocytes. ACS Cent Sci. 2021; 7:1073–88. [PubMed]
  • 30. Wang Q, Tang H, Luo X, Chen J, Zhang X, Li X, Li Y, Chen Y, Xu Y, Han S. Immune-Associated Gene Signatures Serve as a Promising Biomarker of Immunotherapeutic Prognosis for Renal Clear Cell Carcinoma. Front Immunol. 2022; 13:890150. [PubMed]
  • 31. Zhang X, Liu X, Xiong R, An HX. Identification and validation of ubiquitin-proteasome system related genes as a prognostic signature for papillary renal cell carcinoma. Aging (Albany NY). 2022; 14:9599–616. [PubMed]
  • 32. Ribback S, Cigliano A, Kroeger N, Pilo MG, Terracciano L, Burchardt M, Bannasch P, Calvisi DF, Dombrowski F. PI3K/AKT/mTOR pathway plays a major pathogenetic role in glycogen accumulation and tumor development in renal distal tubules of rats and men. Oncotarget. 2015; 6:13036–48. [PubMed]
  • 33. Xuan Y, Chen W, Liu K, Gao Y, Zuo S, Wang B, Ma X, Zhang X. A Risk Signature with Autophagy-Related Long Noncoding RNAs for Predicting the Prognosis of Clear Cell Renal Cell Carcinoma: Based on the TCGA Database and Bioinformatics. Dis Markers. 2021; 2021:8849977. [PubMed]
  • 34. Zhou X, Yao L, Zhou X, Cong R, Luan J, Wei X, Zhang X, Song N. Pyroptosis-Related lncRNA Prognostic Model for Renal Cancer Contributes to Immunodiagnosis and Immunotherapy. Front Oncol. 2022; 12:837155. [PubMed]
  • 35. Dou Q, Gao S, Gan H, Kang Z, Zhang H, Yang Y, Tong H. A Metastasis-Related lncRNA Signature Correlates With the Prognosis in Clear Cell Renal Cell Carcinoma. Front Oncol. 2021; 11:692535. [PubMed]