The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) database was used to download RNA sequencing (RNA-seq) data and clinical information. Additionally, we downloaded gene transfer format files from Ensembl (http://asia.ensembl.org) to identify the lncRNAs from mRNAs. The limma package (v3.40.2) (45) was used to identify differentially expressed lncRNAs between normal and tumor tissues with p < 0.05 and |log₂ fold change (FC)| > 1. Next, we downloaded a list of immune-related genes (IRGs) from the ImmPort database (http://www.immport.org ), which were used in co-expression analysis to identify immune-related lncRNAs (irlncRNAs). Immunogenic cell death-related lncRNAs (icdrlncRNAs) were identified using co-expression analysis of 33 recognized ICD-related genes (ICDRGs) extracted from a large-scale meta-analysis (46). Immune-related and immunogenic cell death-related lncRNAs (IICDLs) were generated from the intersection of two DElncRNA datasets. Thereafter, Spearman’s correlation analysis was used to describe gene–gene correlation. Gene intersection analysis was performed using the ggplot2 package (v3.3.3).

A univariate Cox regression analysis was used to screen for prognosis-related IICDLs. A cluster analysis was performed using the package ConsensusClusterPlus (v1.54.0) to identify IICDLs related to molecular subtypes. Heat map clustering was performed using the pheatmap package (v1.0.12). Kaplan–Meier (KM) analysis was used to compare the prognosis of both clusters. Stack graphs were used to visualize the correlation between clusters and clinical parameters, and chi-square tests were used to analyze them. Visualization was performed using the ggplot2 package (v3.3.3).

The limma package was used to analyze the differential expression of mRNA between clusters 1 and 2. The threshold for differential mRNA expression was defined as p < 0.05 and |log₂ FC| > 0.2. To analyze the biological functions of potential mRNAs, the ClusterProfiler package (v1.54.0) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. The ggplot2 package (v3.3.3) was used for visualization.

We used the prognostic IICDLs to construct a prognostic model using LASSO Cox regression analysis via the glmnet package (v4.1-2). Based on the regression coefficient (β) derived from the multivariate Cox regression analysis, the following formula was used to construct a prognostic signature: risk score = (β_lncRNA1* expression level of lncRNA₁) + (β_lncRNA2* expression level of lncRNA₂) + … + (β_lncRNAn* expression level of lncRNA_n). Each patient was assigned a risk score using this formula in TCGA-HCC cohorts. Patients were categorized based on the median risk score into low- and high-risk subgroups, and overall survival (OS) times were compared between the two groups using KM analysis. A time-related receiver operating characteristic (ROC) analysis was performed using the timeROC package (v0.3) to evaluate the prognostic ability of the risk model.

Both univariate and multivariate Cox regression analyses were performed to test whether the prognostic models are independent of conventional clinical characteristics. A nomogram was developed based on all independent prognostic factors to assess the 1-, 3-, and 5-year survival of patients with HCC (47). Calibration plots were constructed as part of an internal validation process to ensure the predictive accuracy of the nomogram. Nomogram performance was also evaluated using time-dependent ROC analysis. Decision curve analysis (DCA) was performed to determine the clinical net benefit (48).

Data on the somatic mutations in the HCC samples were obtained from TCGA GDC Data Portal in a mutation annotation format. The Maftools package (v2.6.05) was used to visualize and summarize the mutated genes using waterfall plots.

TIMER (https://cistrome.shinyapps.io/timer/) is a web-based resource for estimating the abundance of tumor-infiltrating immune cells (B cells, CD4 and CD8 T cells, neutrophils, macrophages, and dendritic cells) (49). TIMER deduces the abundance of TIICs by deconvolution, Monte Carlo simulations, orthogonal estimates TIMER deduces the abundance of tumor-infiltrating immune cells from gene expression profiles based on a deconvolution method validated by Monte Carlo simulations, orthogonal estimates from DNA methylation-based inferences, and pathological assessments. The abundance of immune cells in each tumor sample was analyzed using TIMER. Eight common immune checkpoint genes (ICGs; SIGLEC15, TIGIT, CD274, HAVCR2, PDCD1, CTLA4, LAG3, and PDCD1LG2) were compared across molecular subtypes and risk groups. Higher Tumor Immune Dysfunction and Exclusion (TIDE) scores were associated with a shorter survival and poorer immune checkpoint blocking (ICB) treatment response. Using the TIDE database, we calculated the TIDE scores for molecular subtypes and risk groups in TCGA. Microsatellite instability (MSI) may be a predictive biomarker of immunity to immune checkpoint inhibitors (50). Therefore, the MSI score was calculated for each sample in TCGA-HCC cohorts to compare the high- and low-risk patients.

The stemness of cancer cells has recently been recognized as a valuable predictive or prognostic factor (51–53). Thus, we used one-class logistic regression (OCLR) to compare stemness for clusters and risks in TCGA (54). Based on the half-maximal inhibitory concentration (IC₅₀) of commonly used chemotherapeutics, including doxorubicin (55), 5-fluorouracil (56), gemcitabine (57), and sorafenib (58), we investigated the molecular subtypes and risk groups associated with them.

The Cancer Cell Line Encyclopedia (CCLE) (https://portals.broadinstitute.org/ccle) database was used to obtain the cell line IICDL expression matrix of tumors. IICDL expression in the HCC tissues and healthy liver tissues was analyzed in TCGA (https://portal.gdc.cancer.gov/), lnCAR (http://lncar.renlab.org/ ), and Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) databases.

The Student’s t-test was used to compare gene expression between tumor and adjacent non-tumor tissue, and the Mann-Whitney U test was used to evaluate proportional differences. K-M analysis was performed to compare the OS between the risk groups and the subtypes via the log-rank test. A multivariate and univariate Cox regression analysis was used for the independent prognosis analysis. Wilcox test was used to compare immune scores between two groups. Statistical analyses were conducted using the R (v4.0.3) software. A significance level of 0.05 was used for all tests. We did not adjust P values for multiple testing (59, 60).

Figure 1 shows the flow diagram of this study. TCGA database, which contains 50 standard samples and 371 tumor samples, was used to obtain the HCC RNA-seq data. Based on the Ensembl gene annotation file, 16,013 lncRNA RNA-seq were obtained, and 147 differentially expressed lncRNAs (DElncRNAs) were initially identified, 25 of which were upregulated and 122 downregulated ( Figure 2A ). There were 143 differentially expressed irlncRNAs (DEirlncRNAs), 24 of which were upregulated and 119 downregulated ( Figure 2B ). Moreover, of the 81 differentially expressed icdrlncRNAs (DEicdrlncRNAs), 16 were upregulated and 65 were downregulated ( Figure 2C ). The intersection of these two DElncRNA datasets was defined as IICDLs, comprising 81 lncRNAs ( Figure 2D ).

The univariate Cox analysis revealed 20 IICDLs associated with prognosis, with AL606489.1 and AC0794666.1 acting as protective factors and the remaining 18 acting as risk factors ( Figure 3A ). According to correlation analysis, most genes were related to one another ( Figure 3B ). Cluster analysis was performed on the 20 prognosis-related IICDLs. Patients with HCC clustered into two subgroups showed the best cluster effect, with good subgroup internal consistency and stability ( Figures 3C–F ). Heatmaps showed that the two clusters significantly differed in gene expression of IICDLs ( Figure 3G ). Cluster 2 had a better prognosis than Cluster 1 according to the survival analysis ( Figure 3H ). In Cluster 2, prognosis was significantly associated with male sex (p < 0.05), better T stage (p < 0.05), early clinical stage (p < 0.01), and early pathological grade (p < 0.05) ( Figures 4A–F ).

We also analyzed differentially expressed genes in C1 and C2. Compared to C2, C1 had 1757 downregulated genes and 7901 upregulated genes ( Figures 5A, B ). Among 1794 IRGs, C1 had 569 upregulated genes (e.g., S100P), 713 downregulated genes (e.g., AQP9), and 512 unregulated genes (e.g., SLP1) in comparison to C2 ( Figure 5C ). Among 33 ICDRGs, C1 had 20 upregulated genes (e.g., BAX), one downregulated gene (FOXP3), and 12 unregulated genes (e.g., PIK3CA) in comparison to C2 ( Figure 5D ).

KEGG and GO analyses were performed based on differentially expressed IRGs and ICDRGs. KEGG analysis revealed that the differentially expressed ICDRGs were mainly involved in the IL-17 signaling pathway, necroptosis, and antigen processing and presentation ( Figure 5E ), whereas the differentially expressed IRGs were involved primarily in the cytokine–cytokine receptor interaction, chemokine signaling pathway, and T cell receptor signaling pathway ( Figure 5F ). GO analysis suggested that the differentially expressed ICDRGs were mainly involved in T cell activation, endocytic vesicle, and MHC protein binding ( Figure 5E ), whereas the differentially expressed IRGs were primarily involved in leukocyte migration, the external side of the plasma membrane, and cytokine receptor binding ( Figure 5F ).

We used the TIMER algorithm, the only method that considers tissue specificity when estimating immune cell populations (61), to determine whether there was a difference in immune infiltration between the two clusters. We found significant differences in the CD4⁺ T cells (p < 0.001), neutrophils (p < 0.001), macrophages (p < 0.001), B cells (p < 0.001), and myeloid dendritic cells (p < 0.001), suggesting that C2 exhibited stronger immunosuppression than C1 ( Figure 6A ). The presence of stromal and immune cells in tumor tissues is indicated by stromal and immune scores; the interaction of cancer cells and tumor stroma affects cancer development, facilitates metastasis, and evades immune surveillance. Therefore, we further investigated the immune and stroma score differences between C1 and C2, revealing that C1 had a higher immune score and a lower stroma score than C2 ( Figures 6B, C ). Additionally, we used the ggplot2 package (v3.3.3) to analyze the ICGs in C2 and C1 and found that CTLA4, HAVCR2, LAG3, PDCD1, and TIGIT were downregulated ( Figure 6D , p < 0.001). TIDE scores for C1 were significantly higher than those for C2 (p < 0.0001), indicating that C2 might achieve greater clinical benefit with ICBs ( Figure 6E ).

Stemness, a molecular marker associated with stem cells, has emerged as a valuable predictor or prognostic factor (51–53). A higher stemness in C1 was noted ( Figure 6F ), suggesting a poorer prognosis for patients in C1. Notably, tumor stemness considerably contributes to cancer chemotherapy resistance (62). Therefore, we assessed whether C1 and C2 responded differently in recommending chemotherapeutic drugs for the treatment of HCC, such as doxorubicin (55), 5-fluorouracil (56), gemcitabine (57), and sorafenib (58). C1 was associated with lower IC₅₀ levels for doxorubicin (p < 0.0001), 5-fluorouracil (p < 0.0001), gemcitabine (p < 0.0001), and sorafenib (p < 0.0001) ( Figures 6G – J ), suggesting that chemotherapy may have more beneficial on C1 than on C2. Additionally, we analyzed somatic mutations in C1 and C2 and found that C1 had a high frequency of TP53 (33%), TTN (24%), and MUC6 (20%) ( Figure 6K ), whereas C2 had a high frequency of CTNNB1 (31%), TTN (26%), and TP53 (26%) ( Figure 6L ).

We then built a prognostic model based on 20 prognosis-related IICDLs. Four IICDLs were tested and chosen for the prediction model based on the results of the LASSO regression analysis ( Figures 7A, B ). Risk score models were calculated based on the following algorithm: risk score = (−0.2854) * TMEM220-AS1 + (−0.0614) * LINC02362 + (−0.0418) * LINC01554 + (−0.0114) * LINC02499. Patients with HCC were divided into two groups based on their risk scores. Figure 7C shows the distribution of risk scores, survival status, and gene expression for these four genes. The KM curves showed that the OS for high-risk patients was significantly worse than that for low-risk patients (p < 0.0001, HR = 2.215) ( Figure 7D ).

The area under curves (AUCs) of time-dependent ROC curves for 1-, 3-, and 5-year OS were 0.706, 0.662, and 0.672, respectively ( Figure 7E ), indicating a good predictive performance. The AUCs of the IICDL-based prognostic model were 0.75, 0.77, and 0.77 for the 1-, 3-and 5-year survival times, respectively.

Furthermore, we identified the differences in risk scores between subgroups based on different clinical pathological factors. The risk score was significantly associated with age (p < 0.05), sex (p < 0.001), T stage (p < 0.05), clinical stage (p < 0.05), and pathological grade (p < 0.05) ( Figures 7F – L ).

We performed univariate and multivariate Cox regression analyses to determine whether other traditional clinical characteristics affected the prognostic model. The TNM stage (p < 0.001, HR = 2.535) and risk score (p < 0.001, HR = 4.6) were independent prognostic factors for OS ( Figures 8A, B ). A predictive nomogram was built to assist in the accurate prediction of clinical outcomes ( Figure 8C ). The predicted OS outcomes matched the actual observations more closely, according to the calibration plot for the internal validation of the nomogram ( Figure 8D ). Additionally, the predictive accuracy of the nomogram and individual prognostic factors were compared using time-dependent ROC curves. The AUCs of the nomogram at 1-, 3-, and 5-year OS were 0.674, 0.716, and 0.721, respectively, which were better than those of the models with only one independent factor ( Figures 8E–G ). DCA was used to determine the clinical relevance of these models, with the combined model predicting outcomes with the highest accuracy ( Figures 8H–J ).

Previous research has linked different clusters to immune status, ICB response, stemness, and chemotherapy response. Accordingly, we aimed to determine whether there were any significant differences between the high-risk and low-risk groups. The TIMER algorithm was applied to investigate the relationship between the tumor immune microenvironment and the signature. The high-risk group had more extensive infiltration of immune cells, including B cells (p < 0.0001), CD4⁺ T cells (p < 0.0001), neutrophils (p < 0.0001), macrophages (p < 0.0001), and myeloid dendritic cells (p < 0.0001), than that in the low-risk group ( Figure 9A ). We also compared the immune, stroma, and MSI scores between the risk groups and found that high-risk patients had lower immune and stroma scores ( Figures 9B–D ) and expressed high levels of five immune checkpoint inhibitors (CTLA4, HAVCR2, LAG3, PDCD1, and TIGIT) ( Figure 9E ). A higher TIDE score was obtained in the high-risk group ( Figure 9F ). The OCLR algorithm indicated a higher stemness in the high-risk group ( Figure 9G ).

We also investigated whether the high-risk and low-risk groups responded differently to doxorubicin, 5-fluorouracil, gemcitabine, and sorafenib. The high-risk group had lower IC₅₀ levels for doxorubicin (p < 0.001), 5-fluorouracil (p < 0.001), gemcitabine (p < 0.001), and sorafenib (p < 0.01) ( Figures 9H–K ), indicating that chemotherapy may have a greater impact on high-risk patients.

To verify IICDL expression, we collected data from HCC patients and cell lines. TCGA data revealed a significant decrease in the expression levels of TMEM220-AS1, LINC02362, LINC01554, and LINC02499 in HCC samples. ( Figure 10A ). According to the GEPIA database analysis, LINC01554 and TMEM220-AS1 were downregulated in HCC ( Figure 10B ). CCLE data indicated that LINC02499, TMEM220-AS1, and LINC01554 were expressed at low levels in most HCC cell lines ( Figure 10C ). The lnCAR database further confirmed these results. TMEM220-AS1, LINC02362, LINC01554, and LINC02499 were significantly downregulated in HCC tissues compared to that in the healthy liver tissues (p < 0.0001) ( Figures 10D–G ). Finally, the effect of single genes on HCC prognosis was analyzed using multivariate Cox regression analysis. The results suggest that LINC01554 and LINC02499 could be used as independent risk factors for determining HCC patient prognosis ( Figure 10H ).