A novel PANoptosis-related long non-coding RNA index to predict prognosis, immune microenvironment and personalised treatment in hepatocellular carcinoma

Data sources

RNA sequencing (RNA-seq) data for individuals with HCC were obtained from The Cancer Genome Atlas (TCGA)-LIHC cohort (https://portal.gdc.cancer.gov/repository) along with the mutation data and associated clinical information. Inclusion criteria for study individuals were patients who contained both RNA-seq data and information on survival time and survival status. Strawberry Perl version 5.32.1.1 was used to differentiate and extract the mRNA and lncRNA matrices for the TCGA-LIHC cohort. Based on previous studies, 24 PANoptosis-related genes (PRGs) were included in this study [9, 10, 12, 13, 30–36] (Supplementary Table 1). The present work is based on open-source data and is free from ethical approval.

Identification of PANoptosis-related lncRNAs in HCC

The R package ‘limma’ was utilised to extract the mRNA expression matrix of PRGs from individuals in the TCGA-LIHC cohort. Using co-expression analysis, lncRNAs associated with PRG mRNAs were obtained and identified as PANoptosis-related lncRNAs (Pearson coefficient > 0.4, p < 0.001). The packages ‘ggplot2’ and ‘galluvial’ were utilised to construct a Sankey diagram of the network of PANoptosis-associated lncRNAs and genes. Additionally, ‘limma’ was utilised to extract differentially expressed PANoptosis-associated lncRNAs in the TCGA-LIHC cohort (|Log2 fold change| > 1, FDR < 0.05), which were visualised using the package ‘pheatmap’.

Identification of a PANRI in HCC

All individuals in the TCGA-LIHC cohort were further randomised into training and validation cohorts in a 1:1 ratio. The univariate Cox (uni-Cox) regression algorithm identified PANoptosis-related lncRNAs associated with survival (p < 0.05) and the packages ‘survival’ and ‘pheatmap’ were utilised to draw forest and heat maps. The optimal PANoptosis-related lncRNAs for constructing PANRI were selected using the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm. In this process, ‘caret’ was used to solve data training for classification and regression problems, and the ‘glmnet package’ was used for variable selection for PANRI. Furthermore, the packages ‘ggExtra’, ‘tidyverse’ and ‘ggplot2’ were utilised to analyse the correlation between the expression of lncRNAs and PRGs in the PANRI and plot correlation heat maps. The risk score for individuals with HCC was calculated using the following equation:

$$\text{Risk}\text{score}={\displaystyle {\sum }_{i=1}^n\text{Coefficient}(i)\times \text{Expression}(i),}$$

where Coefficient represents the regression coefficient of the lncRNA used to develop the PANRI. Accordingly, all samples were classified into high- and low-risk subgroups according to the median score of the training cohort.

Validation of the PANRI in HCC

To validate the potential of the developed PANRI, the ‘survminer’ and ‘survivor’ were utilised to plot Kaplan–Meier (K-M) curves and assess survival differences between high- and low-risk subgroups in the training, testing and TCGA-LIHC cohorts, respectively. The ‘pheatmap’ was also used to validate the risk index by plotting risk heat maps, risk distribution and survival status maps for the different cohorts. Moreover, the uni-Cox and multivariate Cox (multi-Cox) algorithms were used to determine the ability of the index to independently predict the clinical outcomes of individuals with HCC. The receiver operating characteristic (ROC) curves were utilised to assess the prognostic predictive value of the index, which was achieved using the ‘timeROC’, ‘survminer’ and ‘survival’ packages. Finally, the differences in the expression of the 24 PRGs in the two risk subgroups were further analysed.

Establishment and evaluation of a nomogram in HCC based on the PANRI

To better evaluate the clinical outcomes of individuals with HCC, we combined the index with clinicopathological features to construct nomograms according to the outcomes of multi-Cox to predict the survival rates of individuals at 1-, 3- and 5-year intervals. Furthermore, the Hosmer–Lemeshow test calibration curves were utilised to assess the predictability and reliability of the nomogram using the ‘rms’, ‘regplot’ and ‘survivor’.

Validation of the PANRI in clinicopathological parameters

To determine the applicability of the constructed risk index to patients with HCC having different clinicopathological parameters, the ‘survminer’ and ‘survival’ packages were utilised to map the K-M curves of the two risk groups in different clinical subgroups (age, gender, tumour grade and stage). Additionally, ‘ComplexHeatmap’ was utilised to visualise the status of different clinical parameters in the two risk subgroups.

Functional enrichment analysis of the PANRI

Gene ontology (GO) involves the analysis of three aspects: cellular composition, molecular function and biological processes. It is one of the most widely used systems for gene annotation [37]. Differentially expressed genes (DEGs) between the two risk subgroups were obtained using the ‘limma’, followed by GO analysis on the DEGs using the ‘clusterProfiler’, ‘org.Hs.eg.db’ and ‘DOSE’ packages. Enrichment results were mapped using the ‘ComplexHeatmap’, ‘ggpubr’, ‘circlize’, ‘RColorBrewer’ and ‘ggplot2’ packages to analyse the enrichment of DEGs between the risk subgroups in terms of cellular components (CC), molecular functions (MF) and biological processes (BP).

Furthermore, Gene Set Variation Analysis (GSVA) was mainly utilised to assess genomic enrichment in the transcriptome [38]. In this study, GSVA was used to obtain the Kyoto Gene and Genome Encyclopedia (KEGG) pathways enrichment in the two risk groups and further analyse the association between the KEGG pathways and lncRNA expressions in the PANRI. The process was implemented using the ‘limma’, ‘GSEABase’, ‘reshape2’, ‘GSVA’, ‘ggplot2’ and ‘pheatmap’.

Association between the PANRI and TMB in HCC

The downloaded simple nucleotide variation data from the TCGA-LIHC cohort were collated using Perl scripts to obtain a TMB data matrix. The ‘limma’ and ‘ggpubr’ packages were utilised to analyse and visualise TMB differences between the risk subgroups. Moreover, K-M curves were utilised to explore survival differences between different TMB subgroups and different risk subgroup combinations. The 20 genes with the highest mutation frequency in the two risk subgroups were then extracted and a corresponding waterfall was plotted using the ‘maftools’ package.

Predictive value of the PANRI in immune characteristics of HCC

To analyse the association between risk scores and immune characteristics, we performed a single sample gene set enrichment analysis (ssGSEA) using the ‘GSVA’ and ‘GSEABase’ packages, which determined the status of immune cell infiltration, calculated the proportion of different immune cells types in the tumour tissue in each sample and obtained the corresponding immune cell score and function score. Furthermore, the ‘pheatmap’, ‘ggpubr ‘and ‘reshape2’ packages were utilised to plot box plots of the differences in scores between different risk groups.

TIMER is a platform for the systematic investigation of immune infiltration in different types of cancer [39]. Using the TIMER2.0 platform, we acquired a matrix of immune cell infiltration data based on six algorithms for the TCGA-LIHC cohort. Additionally, the ‘ggplot2’, ‘ggpubr’, ‘ggtext’, ‘tidyverse’ and ‘scales’ were utilised to analyse and plot correlation bubble plots, thereby determining the correlation between risk scores and immune cells.

The ESTIMATE algorithm can estimate the number of normal and tumour cells in samples [40]. In the present study, the ‘ESTIMATE’ package was utilised to input the number of stromal cells and immune infiltrating cells in each sample in the TCGA-LIHC cohort and obtain the corresponding scores. The sum of the stromal score and the immune score was considered the ESTIMATE score. Finally, the ‘ggpubr’ was utilised to visualise the differences in scores between the two risk groups.

As a key regulatory factor of the immune system, the activation of immune checkpoints can induce tumour-related immune escape [41]. Therefore, we analysed the differences in different immune checkpoints between the two risk subgroups and visualised the outcomes using ‘ggplot2’ and ‘ggpubr’.

The value of PANRI in predicting response to drug therapy for HCC

To investigate the potential value of the risk index in the clinical pharmacotherapy of HCC, the ‘pRRophetic’ package was utilised to calculate the IC50 of some therapeutic drugs in the different risk groups [42]. Furthermore, for drugs with significantly different IC50s in the two subgroups, box plots were plotted using the ‘ggpubr’.

Consensus clustering analysis

The ‘ConsensusClusterPlus’ package was utilised to cluster patients with HCC based on the PANRI, and principal component analysis (PCA) was performed utilising the R packages ‘ggplot2’ and ‘Rtsne’. The prognostic and tumour immune microenvironment characteristics of patients with different HCC clusters were further analysed using K-M curves, ESTIMATE and immune checkpoints. Box plots of differential outcomes were plotted using the ‘ggplot2’, ‘ggpubr’ and ‘reshape2’ packages.

Cell culture

Human HCC cell lines HepG2 were purchased from American Type Culture Collection (ATCC, USA). Cells were grown routinely in DMEM medium (Gibco, USA) containing 10% FBS (Gibco, USA) in a 37°C humidified 5% CO₂ incubator.

RNA extraction and RT-PCR

Total RNA was extracted from cells using DNAaway RNA Mini-Preps Kit (Sangon Biotech, China), and the NanoDrop 1000 was utilised to assess the quality of the extracted RNA. After reverse transcription, a PCR reaction system was prepared, and qPCR analysis was conducted by real-time detection system through SYBR green I dye (Applied Biosystems, USA) detection. The expression levels of AL049840.4 were quantified and calculated with the method of 2^−ΔΔCt. The primer sequence is shown in Supplementary Table 2.

Transfection of small interfering RNA (siRNA)

For the knockdown of AL049840.4, siRNA was transfected into HepG2 cells using Lipofectamine™ 3000 Transfection reagent (Thermo Fisher Scientific, USA). The final concentration of siRNA was 10 nM, and the volume of transfection reagent was used according to the supplier’s manual. The sequence of siRNA is shown in Supplementary Table 2.

CCK-8 assay

The constructed siRNA-negative control (NC) cells and si-AL049840.4 cells were respectively seeded in 96-well culture plates (5 × 10³/well) and placed at 37°C, 5% CO₂ incubator overnight. The CCK-8 Assay Reagent (Bimake, USA) was added to each well, and the assay was conducted according to the manufacturer’s instructions. A microplate reader (BioTek Instruments, USA) was used to measure the absorbance values at 450 nm, and to detect changes in cell proliferation ability of each group.

Cell apoptosis assay

At 48 h after transfection, cells were harvested through trypsinization, and washed twice with PBS (KeyGen BioTECH, China, pH 7.2). The cells were centrifuged at 1000 g/min for 5 min, then the supernatant was discarded and the pellet was resuspended in 1 × binding buffer (Biolegend, USA) at a density of 1.0 × 10⁵ cells per ml. One hundred μL of the sample solution was transferred to a 5 mL culture tube, and incubated with 5 μl of FITC-conjugated annexin V (Biolegend) and 5 μl of PI (Biolegend) for 15 min. 400 μl of 1 × binding buffer was added to each sample tube, and fluorescence of stained cells was measured on FACSCanto II (BD Biosciences, USA) and analyzed with Flowjo.

Transwell assay

The constructed NC, si-AL049840.4 cells were inoculated into the wells of Transwell plate containing serum-free DMEM medium (4 × 10⁴ cells/well). The lower well contains 800 μl of complete medium (DMEM and 10% FBS). After incubation at 37°C for 36 h, remove cells that have not migrated through the well with a cotton swab. Cells in the lower chamber were fixed with 4% methanal for 15 min and stained with 1% crystal violet in 2% ethanol for 10 min, photographed and counted.

Statistical analysis

All statistical analyses were performed using R software (version 4.1.2) and the corresponding R packages. K-M method was utilised to plot the survival curves of different subgroups. The correlation between different continuous variables was assessed by Pearson correlation test. The Wilcoxon test was utilised for comparing two groups. P < 0.05 was considered as statistically significant for a difference.

Data availability statement

The data sets used and/or analysed during the current study are available from the corresponding author.

PANoptosis-related lncRNAs in HCC

Flowchart reflects the ideas and details of this research (Figure 1). Co-expression analysis identified 1476 PANoptosis-associated lncRNAs and their relationship with PRGs was presented using Sankey plots (Figure 2A). The volcano map showed that 1199 PANoptosis-related lncRNAs were differentially expressed in normal and tumour tissues (|Log2 FC| > 1, FDR < 0.05) (Figure 2B), of which seven lncRNAs were downregulated in tumours and the rest upregulated. Finally, the heat map revealed the expression status of the differentially expressed lncRNAs (Figure 2C).

Figure 1. Flowchart of the present research.

Figure 2. PANoptosis-related lncRNAs in hepatocellular carcinoma. (A) Sankey plots of the correlation between PANoptosis-related lncRNAs and PANoptosis-related genes. (B) Volcano plot showing 7 down-regulated and 1192 up-regulated expressed lncRNAs. (C) Heat map showing PANoptosis-related lncRNAs expressed in normal and tumour tissues.

Establishment and validation of PANRI in HCC

All samples were further randomised into training and validation cohorts (Table 1). Using the uni-Cox regression algorithm, we obtained 282 PANoptosis-associated lncRNAs that were significantly associated with survival (p < 0.05), and the prognostic forest plots of these lncRNAs are presented in Supplementary Figure 1. Furthermore, using LASSO regression analysis (Figure 3A, 3B), seven lncRNAs were screened for use in constructing the risk index (Table 2). The risk score for each individual was calculated using the formula: Risk score = MKLN1-AS × (0.663709706) + ELFN1-AS1 × (0.335059408) + AC015871.6 × (−1.131739203) + AL049840.4 × (0.615830687) + AC026369.2 × (0.514830865) + FOXD2-AS1 × (0.288672846) + LINC00501 × (0.608378767). Figure 3C, 3D show the forest plot and expression heat map of the seven lncRNAs used to construct the index. Additionally, the expression correlation of the seven lncRNAs with 24 PRGs is also presented using a correlation heat map (Figure 3E).

Figure 3. Development of a PANRI in hepatocellular carcinoma. (A, B) The LASSO coefficient and partial likelihood deviance of the PANRI. (C) A risk forest plot of the seven lncRNAs used to construct the PANRI. (D) Expression heat map of 7 lncRNAs used to construct PANRI. (E) Heat map of correlations between the expression of PANoptosis-related genes and the seven lncRNAs used to construct the PANRI.

K-M curves indicated that the expression of seven index-related lncRNAs was correlated with the prognosis of HCC (Supplementary Figure 2). We further validated the stability of the index in the training, validation and TCGA-LIHC cohorts. First, we assessed the K-M curves, risk distribution and survival status of the patients in the three cohorts (Figure 4A–4I). The results revealed that the populations in the low-risk subgroup in all cohorts had a significantly better prognosis. Additionally, expression heat maps showed that AC015871.6 was lowly expressed in the high-risk subgroup in the training, validation and TCGA cohorts, while the remaining six lncRNAs were highly expressed in the high-risk subgroup (Figure 4J–4L).

Figure 4. Validation of the PANRI in hepatocellular carcinoma. (A–C) Kaplan–Meier curves for overall survival in the training (N = 185), validation (N = 185) and entire (N = 370) cohorts. (D–F) Risk score distribution in the three cohorts. (G–I) Survival status in the three cohorts. (J–L) Heatmap of the expression of the seven PANoptosis-related lncRNAs in the three cohorts.

Evaluation of the PANRI in HCC

To further assess the prognostic predictive value of PANRI, we first performed multi- and uni-Cox regression analyses. The outcomes indicated that the PANRI was an independent prognostic element with hazard ratios of 1.040 and 1.052 (Figure 5A, 5B). Additionally, the ROC curves revealed that the area under the curve (AUC) for PANRI at 1, 3 and 5 years was 0.801, 0.726 and 0.700 (Figure 5C). Moreover, the ROC curves also showed that PANRI had better predictive efficacy than clinicopathological parameters such as age, sex, tumour grade and tumour stage (Figure 5D–5F). Notably, differential expression analysis showed that 23 of the 24 PRGs were more highly expressed in the high-risk subgroup (Supplementary Figure 3), suggesting that the high-risk population may have a higher level of PANoptosis.

Figure 5. Assessment of the PANRI in hepatocellular carcinoma. (A, B) Risk forest plots for multivariate and univariate Cox regression. (C) ROC curves of the 1-, 3- and 5-year survival in the TCGA cohort. (D–F) Comparison of risk score ROC curves with clinicopathological parameter ROC curves. (G) Tumour stage and risk status were used to construct a nomogram for predicting patient survival. (H) Calibration curves for the nomogram. ^*p < 0.05, ^**p < 0.01 and ^***p < 0.001.

Nomogram for patients with HCC

Based on multi-cox analyses, the tumour stage and risk status were incorporated into the construction of nomograms for the convenient clinical prediction of clinical outcomes for each patient with HCC (Figure 5G). The corresponding scores for each indicator are presented in the nomogram and all scores are summed as a tool to predict patient prognosis. The nomogram estimated the 1-, 3-, and 5-year OS rates for an individual (high-risk, stage I) to be 0.832, 0.643, and 0.484, respectively. Furthermore, the calibration curves of survival indicated excellent agreement between the actual and nomogram-predicted outcomes (Figure 5H).

Association of PANRI with clinicopathological parameters in HCC

As shown in the K-M survival curves (Figure 6A–6H), individuals with HCC of different ages, grades and stages had better survival rates in the low-risk group. Although the K-M curves for women did not differ significantly in the two risk groups, a trend towards a separation of the curves could be seen. The findings demonstrate the broad applicability of PANRI in patients with HCC having different clinical features. Additionally, the heat map of the clinical parameter status of patients in different risk groups revealed that tumour stage and T-stage differed between the risk subgroups (Figure 6I).

Figure 6. Association of the PANRI with clinicopathological features in hepatocellular carcinoma. (A–H) Kaplan–Meier curves stratified by age, gender, tumour grade and tumour stage. (I) Heat map of the distribution of clinical parameters in different risk groups. ^*p < 0.05, ^**p < 0.01 and ^***p < 0.001.

PANRI-related mechanisms

Given the differences in the prognosis of individuals in the different risk subgroups, we performed GO, GSVA and KEGG enrichment exploration to further investigate the possible mechanisms. GO enrichment analysis revealed that in terms of biological process, DEGs in the two risk groups were mainly enriched in the regulation of cell-cell adhesion, positive regulation of cell adhesion and migration of leukocytes. In terms of cellular composition, DEGs were mainly enriched in the external side of the plasma membrane, collagen-containing extracellular matrix and apical part of the cell. In terms of molecular function, DEGs were mainly enriched in the ECM structure, mainly in signalling receptor activator activity and cytokine receptor binding (Figure 7A).

Figure 7. PANRI-based GO and KEGG analysis. (A) GO analysis shows the enrichment of DEGs between the risk groups. (B) Heat map of functional pathway enrichment differences between the risk groups. (C) Heat map of the spearman correlation analysis between the expression of the seven lncRNAs involved in the model construction and tumour-related pathways.

GSVA explored the differences in KEGG pathways between the risk subgroups. The heat map revealed that the functions enriched in the high-risk subgroup included RNA degradation, ubiquitin-mediated protein hydrolysis, RNA polymerase, pyrimidine metabolism, homologous recombination, DNA replication, mismatch repair, nucleotide excision repair and Notch and mTOR pathways. Notably, these functions are widely involved in the evolution of tumour biological behaviour. Additionally, nitrogen metabolism, phenylalanine metabolism, fatty acid generation, primary bile acid biosynthesis, steroid hormone biosynthesis and linoleic acid metabolism were enriched in the low-risk subgroup (Figure 7B). Finally, spearman correlation analysis heat map showed a broad correlation between the expression of these seven PANoptosis-related lncRNAs and signalling pathways (Figure 7C).

Association of the PANRI with TMB in HCC

The mutation waterfall plot shows the top 20 genes with the highest mutation frequencies in the two risk subgroups (Figure 8A, 8B). In the high-risk population, the top five genes with the highest mutation frequencies were TP53 (35%), CTNNB1 (24%), TTN (22%), MUC16 (17%) and PCLO (11%), while in the low-risk population included CTNNB1 (27%), TTN (25%), TP53 (17%), MUC16 (15%) and PCLO (12%).

Figure 8. Association of the PANRI with TMB in hepatocellular carcinoma. (A) Mutation waterfall map showing the 20 genes with the highest mutation frequency in the high-risk group. (B) Mutation waterfall map showing the 20 genes with the highest mutation frequency in the low-risk group. (C) Comparison of TMB between the two risk groups. (D) Kaplan–Meier curves for the TMB subgroups combined with the risk subgroups.

Tumours with high levels of TMB tend to have higher neoantigen levels, which stimulate the proliferation of anti-tumour effector immune cells and serve as a predictive biomarker for the efficacy of immunotherapeutic response to ICIs in some tumours [43]. Although there was no significant difference in TMB levels between the two risk subgroups (Figure 8C), significant differences in survival between the different four subgroups were observed (Figure 8D). These findings suggest that the combination of the PANRI and TMB levels can better predict the prognosis in HCC.

Role of PANRI in predicting immune characteristics in HCC

The TIME has been closely linked to tumour evolution and influences the response of patients with tumours to immunotherapy. Moreover, there is growing evidence that PCD plays an important role in regulating TIME [33]. To further explore the predictive potential of PANRI in TIME, we utilised the ssGSEA algorithm. Box plots indicated that most immune-related functions including antigen-presenting cell co-inhibition and co-stimulation, immune checkpoint, human leukocyte antigen, inflammation promotion, T-cell co-stimulation and T-cell co-inhibition were stronger in the high-risk population (Figure 9A). In terms of immune cells, immature dendritic cells, macrophages, activated dendritic cells, follicular helper T cells, helper T cells, tumour infiltrating lymphocytes and regulatory T cells were highly expressed in the high-risk subgroup (Figure 9B). Bubble plots showed a positive correlation between most immune infiltrating cells and risk scores (Figure 9C). Furthermore, ESTIMATE analysis revealed that while stromal and ESTIMATE scores did not differ significantly between the two risk groups, immune scores were significantly higher in the high-risk population (Figure 9D–9F).

Figure 9. Correlation of the PANRI with TIME in hepatocellular carcinoma. (A) Box plots of differences in immune-related functions between the high- and low-risk subgroups. (B) Box plots of differences in immune cell scores between the high- and low-risk subgroups. (C) Bubble plots illustrate the correlation between immune cells and risk scores. (D–F) Box plots of the differences in the immune cell, stromal cell scores and ESTIMATE scores in the different risk cohorts. (G) Heat map of the differences in immune checkpoints between the two risk subgroups. ^*p < 0.05, ^**p < 0.01 and ^***p < 0.001.

The key targets in the clinical use of ICIs are currently PD-1, PD-L1 and CTLA-4 [44]. To verify the predictive value of PANRI in the treatment response of ICIs, we analysed the association between immune checkpoints and PANRI. The results showed that most immune checkpoints, including PD-1, PD-L1 and CTLA-4, were more highly expressed in the high-risk group (Figure 9G).

Relationship between the PANRI and HCC therapy

To analyse the therapeutic significance of PANRI in clinical drug selection, a sensitivity analysis of different drugs was performed using the ‘pRRophetic’ algorithm. The IC50 values of specific clinical therapeutics differed between the two risk subgroups (p < 0.001) (Figure 10A–10O). Among them, the targeted drugs, namely ruxolitinib, lisitinib, imatinib, dasatinib, tipifarnib, sunitinib and thapsigargin, and the chemotherapeutic drugs, namely gemcitabine, paclitaxel, etoposide, doxorubicin, 5-fluorouracil and vinorelbine, had lower IC50s in the high-risk subgroup. However, contrasting results were observed for erlotinib, a drug targeting the epidermal growth factor receptor (EGFR).

Figure 10. PANRI-based drug sensitivity analysis. (A–O) Box plots show that the IC50 for certain clinical therapeutics differed significantly between the two risk subgroups (p < 0.001).

Cluster typing based on the PANRI

According to the delta area plots and tracing plots of the consensus clustering algorithm, the lowest intra-cluster difference is found when k = 3 (Supplementary Figure 4). Herein, consensus clustering analysis based on the predictive index was used to classify patients into three clusters (Figure 11A). The tSNE and PCA analysis distinguished the distributional profiles of the three clusters (Figure 11B, 11C). Sankey plots showed that the majority of patients in clusters 1 and 2 were in the high-risk subgroup, while those in cluster 3 were in the low-risk subgroup (Figure 11D). K-M curves indicated that patients in cluster 3 had the best prognosis, followed by cluster 1; however, those in cluster 2 had the worst prognosis (p = 0.002) (Figure 11E). Furthermore, ESTIMATE analysis revealed that samples in cluster 1 had higher immune scores than clusters 2 and 3 (Figure 11F). In terms of immune checkpoints, most immune checkpoints were significantly more highly expressed in cluster 1 (Figure 11G), which, combined with the immune score results, suggests that cluster 1 could be a superior population for ICIs treatment. Thus, the consensus clustering analysis based on PANRI may not only determine prognosis but also aid in TIME characteristic identification.

Figure 11. Hepatocellular carcinoma classification based on the PANRI. (A) Patients were divided into three clusters based on the consensus clustering matrix. (B–C) tSNE and PCA analyses of the three clusters. (D) Sankey plot for the risk cohorts and three clusters. (E) Kaplan–Meier curves of the three clusters. (F) Box plots of the differences in immune cell scores in the three clusters. (G) Expression of immune checkpoints in the three clusters. ^*p < 0.05, ^**p < 0.01 and ^***p < 0.001.

Experimental validation in vitro

In PANRI, AL049840.4 belongs to the risk genes and is significantly associated with the expression of 20 PANoptosis-related genes. Additionally, GSVA results showed that the expression of AL049840.4 in HCC was significantly associated with several tumor-related signaling pathways, however, its role in HCC remains unreported. To know the effect of AL049840.4 on HCC, we further validated it in in vitro experiments. RT-PCR results showed that the transfection of siRNA successfully interfered with the mRNA expression of AL049840.4 in HepG2 cells (Figure 12A). Results of CCK-8 assay showed that knockdown of AL049840.4 suppressed the viability of HepG2 cells (Figure 12B). Additionally, results of cell apoptosis assay revealed that knockdown of AL049840.4 promoted apoptosis in HepG2 cells (Figure 12C, 12D). Furthermore, transwell assay showed a statistically significant decrease in the HepG2 cells migration after AL049840.4 knockdown (Figure 12E, 12F). Together, these results suggest that the expression of AL049840.4 is associated with the viability, migration and apoptosis of HCC cells.

Figure 12. Unfavorable impact of AL049840.4 on HCC in vitro. (A) Validation of knockdown efficiency of AL049840.4 expression in HepG2 cells by qRT-PCR. (B) Cell viability of HepG2 cells after silencing AL049840.4 was detected by CCK-8 assay. (C, D) Cell apoptosis assay revealed that knockdown of AL049840.4 promoted apoptosis in HepG2 cells. (E, F) Transwell assay showed a statistically significant decrease in the HepG2 cells migration after AL049840.4 knockdown.