The workflow of this study is shown in Figure 1 . The transcriptome data of 424 samples (50 non-cancerous samples and 374 tumor samples), the clinical data of 377 samples, and the mutation data of 368 samples of patients with liver hepatocellular carcinoma (LIHC) were downloaded from The Cancer Genome Atlas (TCGA) database. Additionally, the transcriptome data of 445 samples (202 non-cancerous samples and 243 tumor samples) and the clinical data of 260 samples in the International Cancer Genome Consortium-Liver cancer-Riken-Japan (ICGC-LIRI-JP) cohort were obtained from the ICGC database. Healthy samples and any samples with an unclear survival status or a survival time of < 30 days were excluded. After matching the transcriptome data from both databases with eligible survival data, 343 samples from TCGA-LIHC cohort were used for modeling and internal validation. Meanwhile, 230 samples of the ICGC-LIRI-JP cohort were used for external validation. Next, this study obtained a PANoptosis-related gene set comprising 25 pyroptosis-related genes, 8 necroptosis-related genes, and 32 apoptosis-related genes ( Table S1 ) from previous studies (19). To objectively evaluate the differential expression levels of PANoptosis-related genes between liver cancer tissue and non-cancerous liver tissues, the transcriptome data of 110 normal liver tissue were extracted from the University of California Santa Cruz Xena’s Genotype-Tissue Expression(GTEx) project to narrow the gap in sample size between the two tissues.

The PANoptosis molecular subtypes were identified using the consensus clustering algorithm based on the expression of the PANoptosis-related genes in TCGA-LIHC datasets. The clustering process was set to 50 iterations, and 80% of the sample data in each iteration were subsampled to identify a stable and reliable typing. To assess the validity of the classification, visualization was performed using principal component analysis (PCA). The single-sample gene set enrichment analysis (ssGSEA) was used to quantify the abundance of immune cell infiltration (on a scale of 0 to 1). Additionally, the levels of infiltrating immune cells in the three molecular subtypes were compared using the Kruskal-Wallis test. Gene set variation analysis (GSVA) was used to assess the biological functions associated with the molecular subtypes of PANoptosis. The screening criteria to identify the significantly enriched pathway between each two subtypes were as follows: |log fold-change (logFC)| > 0.1 and P (false discovery rate (FDR)) < 0.05. Differentially expressed genes (DEGs) between the molecular subtypes of PANoptosis were identified based on the following criteria: [logFC| ≥ 0.585 and FDR < 0.05.To understand the potential function, DEGs were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene ontology (GO) enrichment analyses (the criterion for significant enrichment was q value < 0.05).

Univariate Cox regression analyses were performed with the expression levels of DEGs in 343 samples of TCGA-LIHC cohort as independent variables and survival time and status as dependent variables (p < 0.05 indicates that DEGs are associated with prognosis). TCGA-LIHC cohort combined with prognostic-related DEGs was randomly split into training and test cohorts in the ratio 1:1. The least absolute shrinkage and selection operator (LASSO) algorithm (when partial likelihood deviance is at its lowest) and multivariate Cox regression analysis were used to further filter prognosis-related DEGs in the training cohort. Finally, a PANRS was obtained. The LASSO method is a compression estimator that can generate finer models by constructing penalty functions. The advantages of the LASSO method include compressing the number of coefficients and reducing multicollinearity problems in regression analysis (20). This study developed a new risk score for patients with HCC based on a gene signature. The risk score was calculated as the sum of the weighted expression of individual genes as follows: risk score = Σ (expression (Gene_n) × coefficient (Gene_n)) (where expression (Gene_n) represents the expression of a specific gene and coefficient (Gene_n) is its corresponding coefficient). The median risk score of the training cohort was used as the optimal threshold to classify TCGA-LIHC, training, test, and ICGC-LIRI-JP cohorts into high-risk and low-risk groups. To evaluate the prognostic accuracy of the signature, receiver operating characteristic (ROC) curve analysis was performed, and the C-index was calculated. The potential clinical relevance of this signature was determined using a clinical applicability analysis. Additionally, a nomogram was developed to further evaluate the accuracy of the signature using a calibration curve.

The cell-type identification by estimating relative subsets of RNA transcripts algorithm was used to obtain the abundance scores of 22 infiltrating immune cells for each TCGA-LIHC sample (total score was 1; p < 0.05 indicates that the accuracy of the levels of infiltrating immune cells obtained in a sample is good). The correlation of the risk scores or risk groups with the abundance of accurate immune cell infiltration was determined using Spearman’s or Wilcoxon tests, respectively. The stromal and immune scores, as well as the combined scores of the two components, in the tumor immune microenvironment were determined using the Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) algorithms. The scores of the high-risk and low-risk groups were compared using the Wilcoxon test. Tumor immune dysfunction and exclusion (TIDE), exclusion, and dysfunction, and microsatellite instability (MSI) scores were used to predict the response of HCC to immunotherapy using the TIDE website (http://tide.dfci.harvard.edu/). Gene expression and clinical information of the immunotherapy cohort are stored in a public database (IMvigor210) (21). The effect of the signature gene lysophosphatidylcholine acyltransferase 1(LPCAT1) on the prognosis of the bladder cancer immunotherapy cohort was evaluated using the Comprehensive Analysis on Multi-Omics of Immunotherapy in Pan-cancer online website (http://camoip.net/). Tumor mutation burden (TMB) can predict immune checkpoint inhibitor (ICI) treatment response as it contributes to the production of neoantigens that activate anti-tumor immune responses (22). Additionally, the neoantigen, as a tumor-specific mutant peptide, is presented only by major histocompatibility complex (MHC) molecules (23). MHC-1 was reported to predict durable ICI response (24). Therefore, this study examined the expression of human leukocyte antigen (HLA) molecules in different risk groups. The variations in each parameter among different groups were analyzed using Spearman’s correlation analysis. The half-maximal inhibitory concentration (IC50) of drugs was predicted using the R package “pRRophetic.”

HCC cell lines (HCCLM3 (BNCC338460), MHCC97-H (BNCC359345), and HepG2 (CL-0103)) and hepatic epithelial cells (THLE-2 (BFN60808733)) were obtained from BeNa Culture Collection (Suzhou, China), Procell (Wuhan, China) and Qingqi Biotechnology Development Co., Ltd (Shanghai, China). These cell lines were cultured in complete medium of Dulbecco’s modified Eagle’s medium. The HepG2 and HCCLM3 cells were transfected in 6-well plates (NEST Biotechnology, Wuxi, China) using Lipofectamine 2000 (Invitrogen), following the manufacturer’s instructions. The corresponding small-interfering RNA (siRNA) sequences are listed in Table S2 . HCC cells were subjected to pyroptosis induction using LPS (1 μg/mL) and ATP (100 μM) as in previous studies (25).

Total RNA was isolated from cells using the Steady Pure Quick RNA extraction kit (Accurate Biotechnology, AG21023). The isolated RNA was reverse-transcribed into complementary DNA using the Evo Moloney-murine leukemia virus reverse transcription premix kit (Accurate Biotechnology, AG11728). qRT-PCR was performed in a LightCycler (Roche, Germany) with 2× Universal SYBR Green Fast qPCR Mix (Abclonal, RK21203). The mRNA expression of LPCAT1 was calculated using the 2^−ΔΔCt method. The primer sequences are shown in Table S2. Western blotting was performed a western blot as described previously (26). The following antibodies were used in western blotting analysis: anti-LPCAT1 antibody (1:500, Proteintech, 16112-1-AP), anti-CASP1 (1:6000, Proteintech, 22915-1-AP), anti-CASP3 (1:500, ABclonal, A2156), anti-CASP7 (1:600, Proteintech, 27155-1-AP), anti-CASP8 (1:500, ABclonal, A0215), anti-GSDME (1:5000, Proteintech, 13075-1-AP), anti-cleaved GSDME antibody(1:500, ABclonal, A23072), anti-GSDMD (Full Length+N terminal, 1:500, ABclonal, A20197), anti-MLKL (1:500, ABclonal, A21894), anti-phospho-MLKL antibody(1:50, ABclonal, A1244), and anti-GAPDH antibody (1:10000, Proteintech, 60004-1-Ig).

Cell proliferation was evaluated using ethynyl-2-deoxyuridine (EdU) staining. Cells from the experimental and control groups were seeded at a density of 20,000 cells per well (100 µL) in a 96-well plate. The EdU-positive cell experiments were conducted using the EdU Apollo 488 kit (RiboBio, C10310-1). EdU‐positive cells rate (%) was calculated as follows: EdU-positive cell rate (%) = [number of EdU‐positive cells (green)/number of Hoechst‐33342-positive cells (blue)] × 100%.

Cell proliferation was examined using the cell counting kit-8 (CCK-8). Cells from the experimental and control groups were seeded in a 96-well plate at a density of 4000 cells (100 µL) per well. The proliferation of HepG2 and HCCLM3 cells at 0, 24, 48, and 72 h post-seeding was determined using CCK-8 reagent (GK10001, GLPBIO, Montclair, USA). The optical density at 450 nm of the reaction mixture was determined using an automatic microplate reader.

The migration ability of HCC cells was assessed using the wound-healing assay. Cells from the experimental group and the control groups were seeded in a 6-well culture plate. When the confluency of cells reached approximately 90%, a scratch was introduced in the monolayer using a 200-µL micropipette tip. The wound closure rate, the following formula was quantified as follows: Wound closure rate (%) = [(original wound area − unhealed wound area)/original wound area] × 100%.

Cell migration was examined using the Transwell assay. HCC cells were trypsinized, and the cell suspension was centrifuged at 1000 rpm for 5 min. The cell pellet was resuspended in phosphate-buffered saline and centrifuged. Next, the cell pellet was resuspended in a serum-free medium. The appropriate amount of cell suspension (300 μL) was inoculated into the Transwell chamber (20,000 cells per well for HCCLM3 and 100,000 cells per well for HepG2) and cultured for 24 h. Non-migrated cells were removed by gently scraping with a cotton swab. The migrated cells were immobilized by using 4% paraformaldehyde, stained with 0.2% crystal violet, imaged, and counted using ImageJ software.

Statistical analysis was performed using R software v4.1.3 and GraphPad Prism. Means between two groups were compared using the Wilcoxon test, while those between three groups were compared using the Kruskal test. The Kaplan-Meir method was used to generate survival curves, which were compared using the log-rank test. In vitro experiments were performed in triplicate, and the data were analyzed using Student’s t-test. Differences were considered significant at p < 0.05.

The expression levels of 86.15% (56/65) of the PANoptosis-related genes were significantly different between HCC (n = 374) and non-cancerous tissues (n = 160). Of these, 33 PANoptosis-related genes were upregulated in HCC samples (11 pyroptosis-related genes, 2 necroptosis-related genes, and 20 apoptosis-related genes), whereas 23 PANoptosis-related genes were downregulated (12 pyroptosis-related genes, 4 necroptosis-related genes, and 7 apoptosis-related genes; Figure 2A ). Additionally, the expression levels of most PANoptosis-related genes (84.54% (53/65); 21 pyroptosis-related genes, 4 necroptosis-related genes, and 28 apoptosis-related genes) were correlated with the overall survival(OS) in patients with HCC ( Figures 2B–D ). In general, the high-expression group exhibited a poor prognosis. These findings support further analysis of PANoptosis-related genes.

The intersection genes of TCGA-LIHC datasets and PANoptosis-related gene set were extracted and matched to the transcriptome data of TCGA-LIHC tumor samples to obtain the expression levels of 65 PANoptosis-related genes in 371 tumor samples. These 371 samples were divided into 3 subtypes, namely PANoptosisClusters A, B, and C, using consensus clustering ( Figure 3A ). The PCA plot provided an intuitive and clear depiction of the three distinct expression patterns in the PANoptosisClusters ( Figure 3B ). To further explore the effect of molecular subtypes of PANoptosis on the prognosis of HCC, the different subtypes were subjected to survival and immune infiltration analyses. The OS of patients with HCC significantly varied according to the molecular subtypes of PANoptosis. PANoptosisCluster B exhibited the best OS, followed by PANoptosisClusters C and A ( Figure 3C ). Additionally, patients in the PANoptosisCluster A exhibit a tumor microenvironment characterized by increased immune cell infiltration. Consistently, PANoptosisCluster B with the best OS exhibited the lowest infiltration of most immune cells. In addition to tumor-promoting immune cells (i.e. myeloid-derived suppressor cells and mast cells), the infiltrating immune cells included anti-tumor immune cells (CD8 (+) T cells and natural killer (NK) cells) and immune cells with both pro-tumor and anti-tumor effects (CD4 (+) T cells and dendritic cells) ( Figure 3D ). Recent single-cell transcriptome sequencing studies on HCC have confirmed that different immune cells express different PANoptosis-related genes (CRADD is upregulated in tumor-specific infiltrating regulatory T cells; TNF is upregulated in cytotoxic FGFBP2+ double‐positive T cells). The same PANoptosis-related gene can be upregulated in multiple immune cells (e.g. NK cells and GNLY+ T cells (cytotoxic CD4T cells) exhibit upregulated expression of GZMB) (27–29). This suggests that immune cells may express several PANoptosis-related genes and that the levels of infiltrating immune cells can indicate the expression level of PANoptosis-related genes to a certain extent. This can explain the high abundance of immune cell infiltration in PANoptosisCluster A with higher expression levels, while the low abundance of immune cell infiltration in PANoptosisCluster B with lower expression levels.

Next, the biological functions of different PANoptosis molecular subtypes were examined using GSVA. The top 10 enriched pathways between PANoptosisClusters A and B, between PANoptosisClusters A and C, and between PANoptosisClusters B and C are shown in Figures 3E–G , respectively. The enrichment of progesterone-mediated oocyte maturation and oocyte meiosis in PANoptosisCluster A was higher than that in PANoptosisCluster B. Additionally, the enrichment of drug metabolism cytochrome p450 in PANoptosisClusters B and C was higher than that in PANoptosisCluster A. PANoptosisCluster B was enriched in olfactory transduction, while PANoptosisCluster C was enriched in mTOR and other signaling pathways. Progesterone-mediated oocyte maturation and oocyte meiosis, which were significantly enriched in PANoptosisCluster A, are closely correlated to HCC progression (30, 31). However, cytochrome P450, an important biosynthetase in drug metabolism, enriched in PANoptosisClusters B and C can inhibit HCC growth by antagonizing HGF/MET signaling or AKT signaling (32, 33). These findings can explain the poor prognosis of PANoptosisCluster A. Additionally, in the olfactory transduction pathway, which was enriched in PANoptosisCluster B, OR1A2 is reported to suppress the proliferation of human HCC Huh-7 cells upon activation with (–)-citronellal (34). The mTOR signal transduction, which was significantly enriched in PANoptosisCluster C, can accelerate HCC progression upon PRIM1 activation (35). These results support that the OS of PANoptosisCluster B is higher than that of PANoptosisCluster C.

A PANRS was constructed to enable the application of these molecular subtypes for the clinical analysis of patients with HCC.Differential analysis of PANoptosisClusters A, B, and C revealed 1565 DEGs ( Figure 4A ). GO and KEGG functional enrichment analyses ( Figures 4B, C ) revealed that the main cellular component in which DEGs were enriched was the chromosomal region. Aberrations in chromosomal regions often lead to HCC cell proliferation, HCC exacerbation, and multiple drug resistance development (36–38). Additionally, the main biological process in which DEGs were enriched was organelle fission. The enhancement of organelle fission is reported to promote HCC metastasis (39, 40) and limit tumor immune surveillance by NK cells (41). Furthermore, the main molecular function in which DEGs were enriched was actin binding. This implies that DEGs can accelerate the progression of HCC by promoting actin binding (42, 43). The KEGG pathway in which DEGs were enriched was human papillomavirus (HPV) infection. Hepatitis B virus infection is a well-known adverse factor for HCC. However, the simultaneous infection of these cells with HPV 16 upregulates the transcriptional activity in HCC (44).

Next, DEGs were subjected to univariate Cox regression analysis that matched transcriptome data and survival information (n=343) and extracted 1155 DEGs associated with prognosis ( Table S3 ). To facilitate subsequent internal cohort validation, TCGA-LIHC cohort (n = 343) was randomly divided into training and test cohorts before establishing the signature. The survival time, survival status, tumor grade and stage, and other important clinical indicators were not significantly different between the two groups ( Table 1 ), indicating that the test cohort can be used as an internal validation cohort.

As detecting a large number of genes in clinical practice is challenging, prognostic DEGs in the training cohort were subjected to LASSO Cox regression analysis to control the number of genes included in the signature. As shown in Figure 4D , the gene coefficient trajectories revealed that an increase in the penalty coefficient (Log λ) results in fewer genes included in the signature, as evidenced by a high proportion of genes with a coefficient of 0. To obtain a signature with an improved fitting effect, the four genes with the smallest cross-validation error in Figure 4E were selected and marked. These genes were subjected to multivariate Cox regression analysis. An optimal signature related to the molecular subtypes of PANoptosis, namely PANRS, was obtained. The PANRS risk score was calculated as follows: risk score = LPCAT1 expression level × 0.387 + chromobox 2(CBX2) expression level × 0.377 ( Figure 4F ). The median risk score of 0.843 was the cut-off value for classifying the train, test, TCGA-LIHC, and ICGC-LIRI-JP cohorts into high-risk and low-risk groups. The Sankey diagram shows the brief process of predicting the prognosis of HCC after the molecular subtypes of PANoptosis were quantified as PANRS ( Figure 4G ). Additionally, comparative analysis of the distribution of risk scores in PANoptosisClusters revealed that PANoptosisCluster A had the highest risk score, followed by PANoptosisCluster C ( Figure 4H ) and PANoptosisCluster B. These results seem to explain the OS of different PANoptosisClusters. Thus, PANRS can well reflect the differential OS between PANoptosisClusters.

The differential OS between the high-risk and low-risk groups in the training cohort was analyzed. Patients in the high-risk group exhibited markedly decreased OS (p < 0.001, Figure 5A ). The risk curve results indicated that patients in the high-risk group were associated with increased mortality rates and upregulated expression levels of the poor prognostic genes LPCAT1 and CBX2 ( Figure 5B ). In the training cohort, the performance of PANRS to predict the OS of patients with HCC was evaluated using time-dependent ROC curves. The area under the curve (AUC) values for predicting 1-year, 3-year, and 5-year OS were 0.793, 0.720, and 0.709, respectively ( Figure 5C ). Similar results were obtained in the test and TCGA-LIHC cohorts ( Figures 5D, E, G, H ). Additionally, the AUC (time-dependent ROC curve) values of PANRS in both test and TCGA-LIHC cohorts were > 0.7 ( Figures 5F, I ), indicating the high prognostic prediction performance of PANRS. Finally, independent prognostic analyses of risk scores in the three cohorts were performed to exclude the interference of other clinical factors. The risk scores were identified as independent prognostic factors for HCC (p < 0.01, Figures 5J–L ).

To assess the potential clinical utility of PANRS, data of TCGA-LIHC cohort were subjected to an applicability analysis. When factors, such as age, sex (male), tumor grade, and tumor stage were considered, the OS of the high-risk group was significantly worse when compared with that of the low-risk group (all p < 0.01, Figures 6A–H ). This suggests that PANRS can be widely used to predict the survival of HCC patients with HCC. To individually determine the probability of 1-year, 2-year, and 3-year survival rates for each patient, a nomogram was established by combining clinicopathological features and risk scores. As shown in Figure 6I , the comprehensive risk score of a patient was 246 points, and the 1-year, 2-year, and 3-year survival probabilities were predicted to be 62.3%, 37.1%, and 28.7%, respectively. The calibration curve in Figure 6J revealed a high level of consistency between the predicted 1-year, 2-year, and 3-year survival probabilities and actual survival rates. Moreover, the C-index and AUC values of the nomogram for predicting 1-year, 2-year, 3-year, and 5-year OS of patients with HCC were higher than those of age, gender, tumor grade, and tumor stage ( Figures 6K–O ). These results indicate that PANRS is a reliable predictive tool with better predictive performance than common clinical indicators. In particular, the individualized predictive performance suggests the potential clinical value of PANRS.

External validation was performed to clarify the generalizability of PANRS. The ICGC-LIRI cohort was divided into high-risk and low-risk groups based on the median risk score of TCGA-LIHC training cohort. As shown in Figure 7A , the high-risk group exhibited poor OS (p < 0.001) and upregulated expression of poor prognosis-related genes (LPCAT1 and CBX2). The AUC values of PANRS for predicting the 1-year, 2-year, and 3-year OS of patients with HCC were approximately 0.7 ( Figure 7B ). Independent prognostic analysis, clinical applicability analysis, nomogram, and calibration charts confirmed that the risk scores can independently, accurately, and individually predict the survival of patients in the ICGC-LIRI cohort for a wide range of populations (all p < 0.05, Figures 7C–F ). Clinical correlation analysis revealed that the expression of the unfavorable prognostic gene LPCAT1 in patients with stage III–IV tumors was markedly higher than that in patients with stage I–II tumors ( Figure S2A ). Additionally, expression of CBX2 was upregulated in patients with a positive family history of tumors in first-degree relatives (p = 0.027, Figure S2B ).

This study also compared the prognostic performance of PANRS with that of previously reported signatures by determining the AUC values. The AUC values of PANRS for predicting the 1-year, 2-year, 3-year, and 5-year OS were higher than those of amino acid metabolism-related gene signature (45), basement membrane-related gene signature (46), and bile acid-related prognostic signature (47) ( Figures S2C–F ). This indicates that the predictive power of PANRS is higher than that of previously reported signatures.

Based on the prognostic power of the risk score for HCC determined in this study, the prognostic power of PANRS for 31 other tumors in TCGA database was examined. The median risk score of PANRS was used to divide TCGA-KIRP cohort into high-risk (n = 196) and low-risk (n = 82) groups. The OS of the high-risk group was significantly lower OS than that of the low-risk group ( Figure S2G ). The AUC values of PANRS for predicting the 1-year, 2-year, and 3-year OS in patients with KIRP were 0.729, 0.746, and 0.712, respectively ( Figure S2H ). The risk scores of patients with KIRP were combined with their clinical characteristics to obtain a nomogram ( Figure S2I ) for assessing the individual survival probability. The survival probability predicted using the nomogram was highly consistent with the actual survival probability, indicating that this nomogram can accurately predict the survival probability of KIRP patients with KIRP exhibiting different clinical characteristics ( Figure S2J ). Additionally, the C-index and AUC values of the nomogram were higher than those of age, sex, and tumor stage, except in the first year ( Figures S2K–O ). These results suggest that PANRS is a powerful prognostic tool whose application may not be limited to HCC.

The tumor microenvironment (TME) is a multifaceted and dynamic ecosystem that significantly affects cancer progression, as well as the efficacy of both immunotherapy and drug treatments (48). The ESTIMATE algorithm was used to obtain the stromal, immune, and ESTIMATE scores. The immune and ESTIMATE scores of the high-risk group were markedly higher than those of the low-risk group (all p < 0.05, Figure 8A ). This suggests an increased infiltration of immune cells in the TME of the high-risk group. Further analysis revealed that the high-risk group exhibited an increased abundance of memory B cells and M0 macrophages, which were both positively correlated with the risk score ( Figures 8B–D ). Conversely, the abundance of naïve B cells, resting memory CD4 T cells, and monocytes was downregulated in the high-risk group. M1 macrophages and monocytes were negatively correlated with the risk score ( Figures 8E, F ). These findings indicated that the high-risk group is highly susceptible to immunosuppressive microenvironments.

ICI, a promising therapeutic for cancer, exerts growth-inhibitory effects against tumor cells by improving the immune function of patients. Only a small proportion of patients benefit from the approved ICIs. The identification of new immune checkpoints can aid in improving the immunotherapy response of patients (49). In this study, the differential expression of 40 immune checkpoint genes (50, 51) in different risk groups was analyzed. The expression levels of 37 immune checkpoints, except for adenosine A2a receptor (ADORA2A), inducible T cell costimulator ligand (ICOSLG), and tumor necrosis factor superfamily 14(TNFSF14), in the high-risk group were higher than those in the low-risk group ( Figure 8G ). Further validation with the TIDE algorithm demonstrated that the high-risk group exhibited decreased TIDE, dysfunction, and MSI scores but increased exclusion scores (p < 0.05; Figures 8H–K ). This suggests a reduced likelihood of tumor immune escape and a potentially favorable response to immunotherapy in the high-risk group.

These findings were further supported by the analysis of the IMvigor210 cohort. The risk scores of the immunotherapy response group (complete remission/partial remission) were markedly higher than those of the immunotherapy non-response group (stable disease/progressive disease) ( Figure 8L ). Additionally, the OS of the high-risk group was higher than that of the low-risk group after immunotherapy (p = 0.011, Figure 8M ), indicating that patients with high risk scores responded well to immunotherapy. The expression levels of the model genes LPCAT1 and CBX2 in the immunotherapy response group were higher than those in the non-responsive group (all p < 0.01, Figures 8N, O ). The upregulated expression of LPCAT1 could also predict a positive response to immunotherapy in patients with bladder cancer (p = 0.032, Figure 8P ).

Additionally, TMB (52) and MHC molecules (24) can predict immunotherapy response. Neoantigens, which are produced by tumors with high TMB, are often associated with improved immunotherapeutic outcomes. However, the downregulation of MHC molecules may suppress the recognition of neoantigens by T cells and consequently decrease the efficacy of immunotherapy (53). This study investigated TMB and MHC expression in different risk groups and demonstrated that increased TMB was positively correlated with the high-risk group. The expression levels of 95.83% (23/24) of the examined HLA genes were significantly upregulated in the high-risk group (p < 0.05, Figures 8Q, R ). These findings suggest that PANRS represents a novel and effective signature for predicting immunotherapy response.

ICIs have improved the clinical outcomes of patients with HCC, increasing the survival rate of patients who responded to the therapy. However, multi-drug resistance mechanisms are the major limiting factors for the efficacy of ICIs (54). Combination therapy, including chemotherapy, is still the main treatment modality. This study compared the IC50 values to predict the sensitivity of HCC populations with different risks to chemotherapeutic drugs. The high-risk group exhibited increased sensitivity to several drugs, including etoposide, cisplatin, gemcitabine, docetaxel, cyclopamine, paclitaxel, pazopanib, rapamycin, sorafenib, and doxorubicin (p < 0.01, Figures 9A–J ). Meanwhile, the low-risk group was sensitive to axitinib, gefitinib, lapatinib, metformin, and AKT inhibitor VIII (p < 0.01, Figures 9K–O ).

In our PANRS, LPCAT1 has a higher risk coefficient. Consider that the sequencing data we previously analyzed were all at the transcriptome level. We therefore assessed the protein expression level of LPCAT1 using the Clinical Proteomic Tumor Analysis Consortium (CPTAC) database. The notable overexpression of LPCAT1 protein in HCC, as shown in Figure S1 , indicates strong consistency between LPCAT1 mRNA and protein expression. However, whether the effect of LPCAT1 on hepatoma cells is consistent with our results still needs to be verified by further in vitro experiments. By qRT-PCR and western blotting, we found that compared with normal hepatocytes THLE-2, mRNA and protein of LPCAT1 were substantially overexpressed in HCC cells (MHCC97-H, HepG2, and HCCLM3) (all p < 0.05, Figures 10A, B ). Furthermore, RNA interference was performed on HepG2 and HCCLM3 with high expression levels. As shown in Figures 10C, D , we successfully inhibited the expression of LPCAT1 protein in HCC cells. By EdU staining and CCK-8 assay, the results indicated that interference with LPCAT1 expression significantly inhibited the proliferative activity of HCC cells (all p < 0.05, Figures 10E–H ). Further analysis using wound-healing and transwell assays revealed a significant decrease in HCC cell migration following interference with LPCAT1 expression (all p < 0.05, Figures 10I–K ). These findings indicate that LPCAT1 contributes to the progression of HCC.

To investigate the effect of LPCAT1 knockdown on pyroptosis, we performed pyroptosis induction in HCCLM3 cells using LPS and ATP and examined the expression level of cleaved GSDMD. LPCAT1 knockdown in HCCLM3 cells promoted the cleavage of GSDMD, generating bioactive GSDMD-N fragments, which can form membrane pores to initiate cellular pyroptosis ( Figure 10L ). In the canonical inflammasome-mediated pyroptosis process, GSDMD is cleaved by activated CASP1, leading to the release of GSDMD-N fragments (55, 56). Consistent with the generation of GSDMD-N, LPCAT1 knockdown upregulated the levels of cleaved CASP1 (p20) ( Figure 10L ). GSDME, a member of the Gasdermin family, is reported to induce pyroptosis under specific conditions (57). LPCAT1 knockdown promoted the cleavage of GSDME ( Figure 10L ). In addition to pyroptosis, LPCAT1 knockdown promoted apoptosis in HCCLM3 cells as evidenced by the minor cleavage of CASP3 (p17), CASP7 (p20), and CASP8 (p18) ( Figure 10M ). Finally, the effect of LPCAT1 knockdown on necroptosis was examined. LPCAT1 knockdown upregulated MLKL phosphorylation in HCCLM3 cells, indicating the induction of necroptosis ( Figure 10N ). These findings indicate that LPCAT1 knockdown promotes PANoptosis in cells.