Prognostic value and immune landscapes of immunogenic cell death-related lncRNAs in hepatocellular carcinoma Open Access

Hepatocellular carcinoma is the sixth most prevalent primary carcinoma and the third most common cause of cancer-related death globally. Hepatocellular carcinoma (HCC) constitutes 90% of liver cancer cases and is one of the most prevalent malignancies with rapid progression and deterioration and an adverse prognosis, with a 5-year survival rate of <50% [1]. Tumor node metastasis classification (TMN) staging is one of the current methods for assessing the prognosis of HCC [2], but its ability to estimate the longer-range prognosis of patients is diminished because it does not account for factors such as cirrhosis [3].

Immunogenic cell death (ICD), a form of apoptosis mediated by activation of immune cells through specific antigen recognition [4], is closely affiliated with tumor development and extensively used in the treatment of tumors (including HCC) [5,6]. In addition, the combination of long noncoding RNAs (lncRNAs) has shown good results as a novel prognostic signature and clinical guidance for patients with cancer [7,8]. ICD-associated lncRNAs have been used in the prognostic assessment of head and neck squamous cell carcinoma, uveal melanoma, low-grade glioma, stomach adenocarcinoma and osteosarcoma [9–13]. The established prognostic signatures demonstrated high predictive accuracy, indicating the importance of ICD-related lncRNAs in the construction of tumor prognostic signatures. Therefore, a novel combination of biomarkers for ICD is needed for the effective prediction of HCC. In the present study, we associated ICD and HCC with lncRNAs to develop a predictive model that is valid and accurate at an early stage, applicable to all stages of HCC and able to provide guidance for personalized treatment after evaluation.

We obtained transcriptomic gene expression information, single-nucleotide mutation information, survival information, and clinical phenotype data from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/, accessed on October 30, 2022) for patients with HCC, yielding a series of 424 specimens, including 374 tumor tissue specimens and 50 normal tissue specimens. The data were processed as shown in Figure 1. Nine tumor tissue specimens with no survival information were excluded. The training cohort and test cohort were obtained after random grouping, and the patient-specific information is shown in Table 1.

A total of 139 ICD-associated genes were obtained from Genecards (https://www.genecards.org/, accessed on October 30, 2022), according to a score>35’ [14]. Pearson correlation analysis (based on |R|²>0.3, P<0.001) was performed by the corr.test function in R (4.2.1).

To identify differentially expressed genes between tumor specimens and normal specimens, we used the ‘DESeq2’ package in R (4.2.1) to conduct differential analysis of lncRNAs obtained from Pearson correlation analysis (|log₂Fold Change|>1 and P<0.05) [15].

Univariate Cox, LASSO regression, and multivariate Cox analyses were performed sequentially in the training cohort by means of the ‘survival’ package in R (4.2.1) [16]. The following formula was used to calculate the risk score for each segment: risk score = (coefficient1 * gene1 expression) +… + (coefficientN * geneN expression). Patients were divided into subgroups according to risk scores; patients with scores above the median were assigned to the high-risk group, and the remaining patients were assigned to the low-risk group.

For the prognostic signature, we performed survival analysis for the training cohort, test cohort, and entire cohort using the ‘survival’ package in R (4.2.1). Time-dependent ROC curves were plotted by using the ‘rms’ package and time-ROC package [17].

To validate the prognostic signature in relation to other clinical information, we used the ‘survival’ package for independent prognostic analysis and the ‘ggDCA’ package, the ‘rms’ package, and the ’foreign’ package for DCA [18].

In accordance with the different clinical characteristics, we calculated the risk scores in different clinical subgroups and grouped them for survival analysis to visualize the applicability of our risk model.

To yield a more comprehensive and accurate risk prediction model, we combined the clinical information of HCC patients and plotted a nomogram and calibration curves using the ‘survival’ package, ‘regplot’ package, and ‘rms’ package [19].

To better visualize the divergent genetic performance among high- and low-risk patients, we performed PCA of complete patient genes, ICD-associated genes, ICD-lncRNA genes, and risk-ICD-lncRNA genes using the ‘scatterplot3d’ package in R [20].

We performed pathway enrichment by GSEA (4.3.2) using the groupings obtained by risk score; specifically, the datasets used were biocarta, kegg, pid, reactome, and go in MSigDB [21]. The pathway enrichment analysis criteria were false discovery rate (FDR) < 0.25 and q-value < 0.01.

We performed immune infiltration analysis of a full set of HCC patients using the CIBERSORT algorithm in R (4.2.1). By using the ‘GSVA’, ‘GSWABase’, and ‘limma’ packages, we performed immune function analysis in both patient groups. We also performed immune scoring using different algorithms (TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, EPIC) through the timer2.0 website (http://timer.cistrome.org/, accessed on November 4, 2022) [22].

From previous papers, we collected information on immune checkpoints and performed an analysis of diverse checkpoint expression in high- and low-risk patients with the ‘reshape2’ package.

By utilizing the ‘maftools’ package in R (4.2.1), we generated waterfall plots of mutation proportions and types and performed analytical calculations of tumor mutation compliance for each of the two patient subgroups according to risk score [23].

We downloaded information on the effectiveness of chemotherapy from the GDSC website (https://www.cancerrxgene.org/, accessed on November 4, 2022) and used the ‘oncoPredict’ package to compare the IC50 differences among 198 chemotherapeutic drugs in combination with information from HCC patients [24].

Eighty-two surgically resected specimens of HCC patients from the Second Affiliated Hospital of Nanchang University were collected. The excised specimens were placed in storage at −80°C after rapid freezing by liquid nitrogen prior to the experiment. We grouped patients and validated the accuracy of the risk gene signature for determining prognosis. Permission for the investigation was granted by the Ethics Committee of the Second Affiliated Hospital of Nanchang University, and informed consent was sought from all patients.

RNA was reverse transcribed at 45°C for 1 h using the SureScript First-Strand cDNA Synthesis kit (GeneCopoeia, China). BlazeTaq SYBR Green qPCR master mix (GeneCopoeia, China) and Applied Biosystems 7500 Fast Real-Time PCR System (Applied Biosystems) were used for subsequent RT-qPCR analysis. The relevant primers used are shown in Supplementary Table S1. The expression of each gene is expressed as the 2^ΔΔCt value. For each specimen, we conducted three experiments to extract the mean value.

To verify the relevance of ICD-associated genes to HCC, we compared the differential expression of ICD-associated genes between HCC and normal tissues through the Human Protein Atlas database (https://www.proteinatlas.org/, accessed on January 25, 2023).

A total of 7222 ICD-lncRNAs were screened out of 14831 lncRNAs by Pearson correlation analysis. Comparing the gene expression levels in 365 tumor tissue specimens and 50 normal tissue specimens, we detected 1006 DELs among 7222 ICD-lncRNAs, of which 906 lncRNAs were up-regulated and 100 were down-regulated (Figure 2A,B).

After univariate Cox analysis of DELs at P<0.05, 243 ICD-associated lncRNAs were obtained (Supplementary Table S2). LASSO regression was used to further screen the relevant lncRNAs (Figure 2C,D), and 11 genes for stable construct models were obtained. We performed multivariate Cox regression analysis on the 11 lncRNAs and finally identified 7 lncRNAs, LINC01136, RP11.298O21.2, RP11.20J15.2, LINC00501, UBE2E1.AS1, LINC01060, and RP5.828H9.1 (Supplementary Table S3 and Figure 2E).

The multivariate Cox analysis results were combined using the following formula: risk score = (1.576275843 × LINC01136 expression) + (−3.014558831 × RP11.298O21.2 expression) + (0.911419424 × RP11.20J15.2 expression) + (1.50737826 × LINC00501 expression) + (4.267512039 × UBE2E1.AS1 expression) + (4.066672785 × LINC01060 expression) + (3.247792567 × RP5.828H9.1 expression). Using this formula, we assigned risk scores to all HCC patients, and we categorized the patients into high- (n=182) and low-risk groups (n=183). The clinical information of patients in each group is shown in Table 2.

Survival analysis was performed for the training cohort, the testing cohort, and the entire cohort (Figure 3A–C), with P-values of <0.001, 0.029, and <0.001, respectively. We subsequently confirmed the correctness of the prognostic signature by time-ROC (Figure 3D–F), demonstrating well-predicted results of the model (training: 1-AUC = 0.852, 3-AUC = 0.777, 5-AUC = 0.751; testing: 1-AUC = 0.649, 3-AUC = 0.648, 5-AUC = 0.607; entire: 1-AUC = 0.705, 3-AUC = 0.691, 5-AUC = 0.689). Heatmaps were drawn for each subgroup of risk-lncRNAs (Figure 3G–I), which consistently showed that RP11.298O21.2 was highly prevalent in the low-risk group, while LINC01136, RP11.20J15.2, LINC00501, UBE2E1.AS1, LINC01060, and RP5.828H9.1 were highly prevalent in the high-risk group. Heatmaps combining clinical details suggested that differentially expressed genes were only associated with risk score and staging and were independent of age and sex (Supplementary Figure S1A). We generated scatter plots of survival status by risk group for the training, testing, and entire cohorts, which showed that risk group was correlated with survival status (Supplementary Figure S1B–G).

Patient age, sex, clinical stage and risk score were used as four factors in the analysis for independent prognosis (Figure 4A,B). The results of univariate Cox regression analysis (HR = 1.806, 1.550–2.06, P<0.001) and multivariate Cox regression analysis (HR = 1.721, 1.466–2.021, P<0.001) indicated that risk score is an individual prognostic factor (Supplementary Table S4). The results of the correlation between the risk score and the three clinical factors indicated that age and sex were not substantially correlated with the risk score, but clinical stage was correlated (Figure 4C–E).

As shown by the survival analysis of different clinical subgroups, the prognostic signature was applicable in all subgroups, especially for patients age <65, female patients, and stage I∼II and T1∼2 patients (Supplementary Figure S2).

PCA was carried out separately for all genes, ICD-associated genes, ICD-associated lncRNAs, and risk-ICD-lncRNAs. 3dPCA plots showed that the seven prognosis-associated lncRNAs were most strikingly different between the high- and low-risk groups (Supplementary Figure S3).

We plotted ROC curves for the four factors (Figure 5A), and the AUCs were 0.705, 0.540, 0.516, and 0.656 for risk score, age, sex, and clinical stage, respectively. The above results illustrate that the nomogram has good predictive efficacy and provides accurate prognosis evaluation results for HCC patients.

In an attempt to yield a more comprehensive and accurate prognostic score, we integrated the patients’ clinical information with the risk score to derive a nomogram (Figure 5B) to predict the patients’ 1-, 3-, and 5-year survival rates. Decision curve analysis (DCA) further confirmed that the model incorporating the risk score was more worthy of evaluation (Figure 5C). The calibration curve demonstrated the predictive capability of the nomogram (Figure 5D,E), suggesting that the comprehensive nomogram with the risk score (C-index = 0.695) is superior to the nomogram without the risk score (C-index = 0.671). Furthermore, we arbitrarily chose one patient to confirm the feasibility of prediction with the nomogram (Supplementary Figure S4).

In the Biocarta enrichment analysis, the low-risk gene set was functionally associated with endogenous thrombospondin activation, lipid metabolism and toxic nuclear receptors, and complement pathways; the high-risk genes were relevant to cell cycle regulation, BRCA1/2, and ATR susceptibility in cancer (Figure 6A). In the remaining four gene sets, the pathways in the low-risk group were closely affiliated with substance or drug metabolism, while the high-risk group was chiefly connected with mitosis (Figure 6B–E). The outputs of the GSEA for the high- and low-risk groups are exhibited in Supplementary Tables S5 and S6.

In the immunoinfiltration analysis, antigen-presenting cells (APCs), macrophages, follicular regulatory T lymphocytes (Tfhs), T helper 2 cells (Th2), and regulatory T lymphocytes (Tregs) were highly frequent in the high-risk group; B cells and natural killer (NK) cells were highly expressed in the low-risk group (Supplementary Figure S5A), as detailed in Supplementary Figure S6A–J.

In the immune function analysis, APC_co_stimulation, HLA, and MHC_class_I scores were higher in the high-risk group, while Type_II_IFN_Response scores were higher in the low-risk group (Supplementary Figure S5B).

The immune differential analysis by seven different algorithms is shown in the significance bubble plot, which had similar results to those of ssGSEA, indicating that immune cells tend to infiltrate more in high-risk patients (Supplementary Figure S5C).

The violin plot illustrates a salient variance in TMB across high- and low-risk patients with a P-value = 0.004 (Figure 7A). The presence of a higher tumor mutational burden in the high-risk group was further confirmed by the TMB percentage bar graph (Figure 7B). In the intuitive demonstration of SNPs, 87.22% of patients in the low-risk group had mutations (Figure 7C), and 89.2% of patients in the high-risk group had mutations (Figure 7D). The top five mutant genes in the low-risk group were CTNNB1, TTN, TP53, MUC16, and ALB, and the top five mutant genes in the high-risk group were TP53, CTNNB1, TTN, MUC16, and ALB.

Among the 35 immune checkpoints with differences, we detected high expression of immune checkpoints (e.g., PD1 and CTLA4) in the high-risk group (Supplementary Figure S6K), implying that it may be possible that the high-risk group may benefit more from immune-associated therapies against these checkpoints.

By comparing the IC50 of anticancer drugs across high- and low-risk patients, we yielded an aggregate table of 91 drugs that could be used in clinical treatment (Supplementary Table S7). High-risk patients were more sensitive to drugs related to mitosis, the cytoskeleton, protein stability destruction and part of the cell cycle (such as MK-1775, AZD7762, and palbiclib), while low-risk patients were more sensitive to DNA replication-related drugs (e.g., nelarabine, oxaliplatin, fludarabine, and mitoxantrone). The remaining 107 drugs with no difference between the high- and low-risk groups are displayed in Supplementary Table S8, indicating that the model cannot determine the group preference of these drugs. A summary of the effective drugs and their differences is shown in Supplementary Table S9.

We derived the protein expression levels of five of the ICD-associated genes in HCC and normal tissues from the HPA database (Figure 8A).

The expression of seven ICD-associated risk lncRNAs was validated in RT-qPCR experiments in high- and low-risk patients. The results demonstrated that all seven lncRNAs were differentially expressed (Figure 8B), and all but RP11.298O21.2 were highly expressed in the high-risk group.

HCC is a global health concern, with the incidence reaching >90,000 cases by 2025 [25,26]. ICD is a modality that is closely associated with tumorigenesis development by recognizing surface antigens that further induce immune cell collections on tumor cells, thereby promoting necrosis through immune cell infiltration [27]. To address the prediction problem, we established a prognostic signature of LINC01136, RP11.298O21.2, RP11.20J15.2, LINC00501, UBE2E1.AS1, LINC01060, and RP5.828H9.1 in HCC and plotted a nomogram to assess patient prognosis. The results of the survival analysis showed that our signature has good assessment ability, and the calibration curve indicated that the risk score added predictive accuracy to the nomogram.

Previous studies have shown that ICD is closely related to tumor development [27]. Hence, many researchers have experimented with ICD-lncRNAs to build prognostic signatures. In the ICD-lncRNA model of HNSC, patients in the high-risk group presented a noticeably dismal prognosis (P<0.0001 in TCGA, P=0.013 in GSE41613) [9]. The same results were demonstrated in ICD-lncRNA models of UM (P<0.01 in TCGA, P=0.009 in GEO) and GBMLGG (P<0.001 in TCGA, P=0.003 in CGGA) [10,11]. Consistent with the prognostic signature constructed in the present study (P<0.001 in the training cohort, P=0.029 in the testing cohort), the prognosis was worse in the high-risk group. In these studies, we noticed that immune recognition sites and killer immune cells were less abundant in the high-risk group. This led to the assumption that ICD was an important process of self-protection against tumors. Collectively, the ICD-associated prognostic signature was a worthy method for evaluating prognosis.

The progression of HCC is strongly related to the dysregulation of the host's immune mechanisms. Alessandro Granito et al. described that Tregs, with CD4+ CD25+ Foxp3+ as molecular markers, had a suppressive effect on CD8+ T cells and NK cells in the development of HCC, thus reducing the effect of the host immune system on tumor cell death [28]. Consistent with the results of the present study, ssGSEA results revealed a significant elevation of Tregs in the high-risk group, accompanied by a decrease in NK cells. Consequently, we speculated that the detrimental effect of ICD on tumor cells in the high-risk group was suppressed by Tregs so that the effect of the immune system on HCC cell death was reduced, leaving the patients in the high-risk group without the protection of ICD and causing a poor prognosis.

The GSEA results indicated that the malignancy of HCC in the high-risk group was related not only to the TME but also to HCC cells. Lu’s study illustrated that c-MYC and cyclin D1 were both targets of the WNT pathway, and up-regulation of β-catenin produced by the WNT pathway led to up-regulation of both targets [29]. In addition, Fubo Ji et al. found that significant up-regulation of P450 could hinder transcription of the downstream gene c-MYC and prevent HCC [30], consistent with our analysis that P450 was up-regulated in the low-risk group and c-MYC was enriched in the high-risk group. Hence, the combined effect of enhancing P450 expression and inhibiting c-MYC expression can be considered to enhance the prognosis of patients.

The results of the drug sensitivity prediction analysis provided guidance for individualized treatment after evaluation. Immunotherapy and chemotherapeutic agents are two of the commonly used options for HCC [31]. Sagar et al. found that ICD could be mediated through the p38 pathway, leading to the effective death of tumor cells [32]. Zhuo et al. showed that directing apoptosis can effectively increase ICD and thus achieve anti-HCC effects [5]. EGFR signaling drugs can inhibit cell proliferation and enhance ICD by increasing the permeation of CD8+ and NK cells for the treatment of HCC [33]. Furthermore, Bernardo Stefanini et al. found that the combination of TKIs with ICIs is expected to be a mainstream regimen in the future [34]. While reducing immunosuppression through immune targeting and strengthening ICD, this combination also suppressed tumor proliferation by inhibiting tyrosine kinases (including EGFR signaling and RTK signaling) and effectively reduced drug resistance during the treatment process. Our drug sensitivity prediction results confirmed that the high-risk group was more sensitive to TKIs than the low-risk group, and the immune checkpoint analysis showed high expression of some popular ICI targets (e.g., PDCD1 and CTLA4) in the high-risk group, which may provide a new direction for personalized treatment of patients in the high-risk group.

Our study established a novel prognostic signature of ICD in combination with lncRNAs in HCC. In the present study cohort, the predictive performance of TMN staging was inferior to that of the obtained prognostic model. Wu et al. established a prognostic signature for ECM-related lncRNA HCC with AUCs of 0.746, 0.683, and 0.670 for time-ROC curves at 1, 3, and 5 years in their TCGA cohort [35]. In a second investigation, Huang et al. developed a genome instability-derived lncRNA-based prognostic signature for HCC with an AUC of 0.619 in a multi-ROC analysis [36]. In the present study, the TCGA cohort time-ROC AUCs of 0.705, 0.691, and 0.689 were significantly higher than those of the proptosis-related prognostic signature, and the AUC of 0.705 for the entire cohort multi-ROC analysis was also better than that of the inflammation-related prognostic signature. We have for the first time aggregated the pathway mechanisms and antitumor drug sensitivities of different risk groups to visualize the potential mechanisms of predictive models and the value of clinical drug guidance. However, it should be acknowledged that our experiments still have shortcomings. First, we had a limited sample size to download from TCGA, and the number of patients with N1-3 and M1 stages was small and all in the high-risk group, so we were unable to perform subgroup survival analysis for N and M staging. Second, we could not obtain BALC staging of HCC patients directly from the TCGA database, so we were unable to compare its predictive performance with that of the model. Finally, we utilized an online database, and subsequent experiments need to be conducted. New data still need to be collected to study the mechanism of these seven ICD-associated lncRNAs.

We constructed an effective prognostic signature for patients with HCC using seven ICD-lncRNAs, which provided guidance for the prognostic assessment and personalized treatment of patients. In addition, ICD was integrated with WNT pathway-induced immune escape and malignant tumor biological behavior for the analysis of potential mechanisms. To this end, we also provide guidance on clinical agents (e.g., mitogenic inhibitors, JNK, and p38 pathway agents) for the treatment of patients after evaluation. Due to the limitations of experiments and specimens, we need more basic, mechanistic and clinical studies for validation.

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

The authors declare that there are no competing interests associated with the manuscript.

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.

Wanying Chen: Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review & editing. Kexin Shu: Conceptualization, Data curation, Formal analysis, Investigation. Chenxi Cai: Conceptualization, Resources, Data curation, Software, Investigation. Jiatong Ding: Conceptualization, Resources, Data curation, Software, Formal analysis, Investigation. Xin Zhang: Conceptualization, Resources, Data curation, Software, Formal analysis, Investigation. Wenxiong Zhang: Conceptualization, Funding acquisition, Validation, Visualization, Methodology, Writing—original draft, Writing—review & editing. Kang Wang: Resources, Software, Writing—original draft, Project administration, Writing—review & editing.

The current study investigated the publicly available data and was also approved by the ethics committee of the Second Affiliated Hospital of Nanchang University (Nanchang, China). All methods were carried out in accordance with the Declaration of Helsinki.

The written informed consent was obtained from all participants included in the study.

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.

95%H

95% high

95%L

95% low

APC

antigen-presenting cell

AUC

area under curve

C-index

Concordance index

CGGA

Chinese Glioma Genome Atlas

Coef

coefficient

DEL

differential expression of long non-coding ribonucleic acid

FDR

false discovery rate

GEO

Gene Expression Omnibus database

GSEA

gene set enrichment analysis

HCC

hepatocellular carcinoma

HPA database

The Human Protein Atlas database

HR

hazard ratio

IC50

half maximal inhibitory concentration

ICD

immunogenic cell death

ICD-associated gene

immunogenic cell death-related gene

ICD-lncRNA gene

immunogenic cell death long non-coding ribonucleic acid gene

iDAMP

inhibitory damage-associated molecular pattern

IQR

interquartile range

K-M survival analysis

Kaplan–Meier survival analysis

LncRNA

long non-coding ribonucleic acid

M

metastasis

N

node

NES

standardized enrichment score

NK cell

natural killer cell

PCA

principal component analysis

Pr>|z|

Probability>|z score|

Risk-ICD-lncRNA gene

risk immunogenic cell death long non-coding ribonucleic acid gene

ROC

receiver operating character curve

RT-qPCR

reverse transcription quantitative-polymerase chain reaction

Se(coef)

standard error (coefficient)

SNP

single-nucleotide polymorphisms

ssGSEA

single sample gene set enrichment analysis

T

Tumor

TCGA

The Cancer Genome Atlas

Tfh

follicular regulatory T lymphocytes

Th2

helper T 2 cells

Time-ROC

time-dependent receiver operating character curve

TMB

tumor mutation burden

TMN

tumor node metastasis classification

Treg

regulatory T lymphocytes

Z

Z-score, standard deviation