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

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.


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

Data sources and access
The University of California, Santa Cruz (UCSC) Xena database (https://xena.ucsc.edu/, 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 (https://portal.gdc.cancer.gov/repository, 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 (https://www.genecards.org/, 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, ICDrelated 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 https://www.gsea-msigdb.org) [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 (http://timer.cistrome.org/, 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 (http://tide.dfci.harvard.edu/, until November 1st, 2022) was used to retrieve the TIDE scores and related data of KIRC patients [17]. Therefore, we conducted singlesample 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 (https://www.proteinatlas.org/) 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).
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).  (LINC01 192). Clinical information for all patients in the highand 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 3A-3C), 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

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

Enrichment analysis
KEGG enrichment analysis revealed the related functions of variably expressed genes, including ubiquitinmediated 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 AGING  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 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.

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

DISCUSSION
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 metastasisrelated 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.

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

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.

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

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