The TCGA database (https://portal.gdc.cancer.gov) was utilized to procure RNA sequencing and clinical information for 374 HCC samples. Additionally, five single-cell RNA sequencing (scRNA) data of HCC were downloaded from GSE242889 (Li et al., 2024) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi) to visualize PCD scores within immune cells. External validation data was sourced from ICGC-LIRI-JP (https://dcc.icgc.org/projects/LIRI-JP) and GSE10143 (Hoshida et al., 2008). Three analytical platforms were employed in the comprehensive analysis. The Metascape online tool (https://metascape.org/) was utilized for the analysis of functional annotations (Zhou et al., 2019), while TIDE (http://tide.dfci.harvard.edu/) was conducted for the prediction of immunotherapy responses (Jiang et al., 2018). Additionally, Proteomaps (https://proteomaps.net/) was performed to illustrate the composition of protein in each groups (Liebermeister et al., 2014). Finally, the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) was consulted for immunohistochemical (IHC) staining evaluation comparing healthy and HCC samples.

RNA-seq data normalization was performed utilizing the “DESeq2” R package, while background correction and normalization of microarray data were done by the “affy” package. ScRNA-seq data was normalized by utilizing the “NormalizeData” function within the “Seurat” package.

PCD-associated genes were obtained from existing literature sources (see Supplementary Table 1). Utilizing the CIBERSORT tool, a deconvolution algorithm was employed to calculate 22 different immune cell types using bulk-seq data (Chen et al., 2018). In order to maintain consistency and accuracy in our gene expression data, we adhered to standard preprocessing procedures for normalizing RNA-seq and microarray data. The TME score was determined by calculating enrichment scores generated by CIBERSORT.

We conducted univariate Cox regression analysis with bootstrap resampling (1,000 iterations) to identify potential prognostic PCDs related to OS in the TCGA-LIHC dataset. To further refine the selection of prognostic PCDs, we employed the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis using the “glmnet” R package. A bootstrap algorithm (resampling = 1,000) was employed in a multivariate Cox regression analysis to identify the most correlated PCDs with prognosis. To the stability of the results, we fixed the bootstrap coefficient of each included PCD: b o o t s t r a p   c o e f f i c i e n t = c o e f f i c i e n t b o o t s t r a p   s t a n d a r d   d e v i a t i o n

The PCD score was calculated using the formula below: P C D   S C O R E = ∑ i = 1 n b o o t s t r a p   c o e f f i c i e n t   i n c l u d e d   P C D i × e x p r e s s i o n   l e v e l   i n c l u d e d   P C D i

We used the median as a cutoff point for categorizing samples into high- and low-score groups. The TCGA–LIHC cohort investigated survival differences among two PCD score groups by the “survival” package.

The TME score was determined by quantifying the presence of immune cells in 22 subtypes of HCC utilizing the CIBERSORT algorithm. Subsequent survival analysis was conducted for individual patients based on the immune cell infiltration, Prognostic immune cells have been identified as those demonstrating differential overall survival rates among various subgroups. Additionally, the prognostic immune cells’ bootstrap coefficient was computed through multivariate Cox regression analysis (resampling = 1,000). The TME score was defined as: T M E   S C O R E = ∑ i = 1 n b o o t s t r a p   c o e f f i c i e n t   p r o g n o s t i c   i m m u n e   c e l l i × i n f i l t r a t i o n   l e v e l   p r o g n o s t i c   i m m u n e   c e l l i

Patients were categorized into two subgroups, TME low and TME high, based on their median TME score. This stratification was done in order to conduct a survival analysis and evaluate differences in OS between the two subgroups. Subsequently, the PCD–TME classifier was formulated by integrating the PCD and TME scores. Within the study, patients with HCC were segregated into four distinct subgroups according to their PCD and TME scores: PCD high-TME low, PCD low-TME low, PCD high-TME high, and PCD low-TME high. To simplify the PCD-TME classifier for clinical application, PCD high-TME high and PCD low-TME low were amalgamated into a single category labeled as Mixed due to their less divergence. The study proceeded with a survival analysis to examine variations in OS across the three groups. Following this, the “timeROC” and “survivalROC” R packages were applied to evaluate the efficacy of the PCD-TME classifier by determining the area under the curve (AUC) of receiver operating characteristic (ROC) curves at 1-, 3-, and 5-year time points.

We performed survival analysis in the TCGA-LIHC cohort to examine variations in OS among subgroups. The PCD-TME classifier was evaluated for its potential as an autonomous prognostic indicator for hepatocellular carcinoma (HCC) within the TCGA-LIHC through Cox regression analyses. Additionally, the findings of this study were corroborated through validation in the ICGC-LIRI-JP.

Gene set enrichment analyses (GSEA) were carried out to identify pathways linked to low PCD-high TME and high PCD-low TME. Genes with comparable expression profiles were clustered using weighted gene co-expression network analysis (WGCNA) (Costanzo et al., 2016; Giulietti et al., 2016) through an unsupervised analysis approach. Subsequently, a Metascape analysis was conducted to visually represent the enrichment results for the genes that were identified as key modules by WGCNA.

An effective anti-tumor immune response can be initiated through tumor mutational burden (TMB) (Meléndez et al., 2018), leading to prolonged clinical outcomes. TMB levels were compared among subgroups by calculating individual TMB scores for samples in the TCGA-LIHC cohort using established methodologies. Subsequently, hub genes exhibiting differences among high PCD-low TME and low PCD-high TME were identified and analyzed. TMB scores for each tumor were also computed following established protocols. The differential gene expression analysis was utilized by the R package of “limma.” Furthermore, Proteomaps were constructed utilizing a web-based tool. The R package “clusterProfiler” was employed to conduct the analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways across all cohorts (Liebermeister et al., 2014).

In the analysis of scRNA sequencing, cells were filtered based on the criteria of detecting between 300 and 5,000 genes, with mitochondrial genes accounting for less than 30% of the total. Data normalization was performed using the LogNormalize method, and the top 2,000 most variable genes were identified for downstream analysis. Dimensionality reduction was conducted using Principal Component Analysis (PCA), and the first 20 components were used to cluster cells. Batch effects were corrected using the Harmony algorithm. The “inferCNV” package was used to identify malignant cells. Additionally, intercellular communication patterns were elucidated using the “CellChat” package.

Following the quantification of protein concentration, the protein sample underwent electrophoretic separation on 12% SDS−PAGE gels and subsequent transfer onto 0.45-mm PVDF membranes (Millipore). Subsequently, the membranes were subjected to incubation with HTRA2 Antibody (AF1855, Beyotime) at 4°C overnight, followed by HRP-conjugated secondary antibodies (AS014, ABclonal). The primary antibody was diluted at a ratio of 1:1,000, while the secondary antibody was diluted at a ratio of 1:3,000. Finally, images were captured using a ChemiDoc Imaging System (Bio-Rad, United States), and quantitative analysis was performed using ImageJ. Student’s t-test was used for data comparisons between JHH-7 HCC cells and L-O2 normal liver cells.

Total RNA was extracted from JHH-7 HCC cells and L-O2 normal liver cells using the RNA Isolation Kit (R0017M, Beyotime). RNA quantity and concentration were measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific, United States). Reverse transcription of total RNA to cDNA was performed using the BeyoRT™ III M-MLV Reverse Transcriptase Kit (D7176M, Beyotime). Quantitative real-time PCR (qRT-PCR) was then conducted using the Taq Pro Universal SYBR qPCR Master Mix Kit (Vazyme, China). The cycling threshold (Ct) for HTRA2 was recorded, and the relative expression of HTRA2 mRNA was calculated using the 2^−ΔΔCT method, with appropriate controls. The primers used in the qRT-PCR protocol are listed in Supplementary Table 2.

Statistical analyses were performed with various methods, including Cox regression analyses, log-rank tests, Wilcoxon rank-sum tests, Student’s t-tests, and Fisher’s exact tests in R 4.1.1. Multiple testing correction was implemented using the Bonferroni method for comparisons involving multiple groups. Statistical significance was set at p < 0.05 unless specially indicated. * P   <   0.05 ; ** P   <   0.01 ; *** P   <   0.001 .

The schematic illustration of the overall research, as shown in Figure 1, involved the evaluation of 374 HCC samples from TCGA-LIHCC to develop a method for indicating PCD expression. Through an initial screening of difference genes (DEGs) and prognostic genes using univariate Cox analysis and a bootstrapping algorithm, a total of 461 prognostic PCDs were identified (Supplementary Table 3). Figure 2A presents the heatmap illustrating the top 20 prognostic PCDs that are up-regulated in HCC. To ascertain the most robust prognostic genes among the candidates, LASSO regression analysis was utilized to determine their risk prediction contributions, as illustrated in Figures 2B, C. Multivariate Cox analysis identified 23 PCD-related genes that significantly influenced PCD scores, as depicted in Figure 2E. Subsequently, patients in the TCGA-LIHC cohort were stratified into two groups on the basis of their PCD scores. Statistical analysis indicated that patients with lower PCD scores performed more favorable clinical outcomes, whereas those with higher PCD scores experienced poorer outcomes (Figure 2D). Furthermore, GSEA revealed a potential distinction between the PCD subgroups, suggesting that tumors characterized by a high PCD phenotype are more proliferative (Figure 2F) and all the enrichment score results could be assessed by Supplementary Table 4. A PCD score was developed to forecast the prognosis of the patients with HCC in this section, and the underlying function of PCD molecules in HCC.

An immune cell signature was derived from the transcriptomes of 22 immune cells utilizing the CIBERSORT algorithm. Six distinct immune cell types were identified based on their optimal cutoff values, demonstrating their protective roles in OS. These immune cell types consist of resting memory CD4⁺ T cells, CD8⁺ T cells, activated NK cells, naïve B cells, resting mast cells, and M1 macrophages (Figures 3A–F). Figure 3G displays the multivariate Cox analysis of immune cells. It was observed that patients with high TME scores had significantly longer survival compared to those with low TME scores, in contrast to the PCD score, as shown in Figure 3H. Furthermore, Figure 3I presents a correlation analysis demonstrating the association between immune infiltration and PCD expression. The complement system is a group of proteins that promote the removal of microorganisms and damaged cells by antibodies and phagocytes, and it enhances (complements) antibody- and cell-mediated immune mechanisms through a series of cascading reactions, which was enriched in the TME high groups (Figure 3J). The result of the GSEA analysis of TME scores is available in Supplementary Table 5. TME scores and elucidation of the relationship between immune cells and PCD were established in this section.

In light of the aforementioned research, the inquiry arose as to whether a PCD-TME classifier could be developed through the integration of PCD and TME scores. This amalgamation yielded three distinct subgroups: PCD high-TME low, PCD low-TME high and Mixed. A notable disparity in prognosis was observed between the PCD-TME classifier and the TCGA-LIHC cohort, as illustrated in Figure 4A. In contrast to the other two groups, samples in the subgroup characterized by PCD low-TME high exhibit a more favorable prognosis. The time-dependent ROC curves were applied to predict the capability of the PCD-TME classifier. The AUC values were 0.747 for 1 year, 0.757 for 3 years, and 0.766 for 5 years, as depicted in Figure 4B. WGCNA was performed to investigate gene variations among these subgroups in order to elucidate significant differences in survival outcomes. The “pink” and “blue” modules were identified as having the most pronounced variations between PCD low-TME high and PCD high-TME low cohorts, as illustrated in Figures 4C, D. All genes corresponding to the “pink” and “blue” modules were subjected to functional annotation using the Metascape. Enrichment analysis indicated significant enrichment of cell cycle-related pathways in the PCD high-TME low group, whereas tissue development-related pathways were predominantly found in the PCD low-TME high group (Figures 4E, F). We also tested the PCD-TME classifiers by the R package of “fgsea,” which shows a similar result (Supplementary Figure 1). A PCD-TME classifier was developed in this section by integrating PCD and TME scores, and further investigation was conducted to explore the functional distinctions between the two subgroups.

Figures 5A, B illustrate that the PCD-TME classifier performed a statistically significant association with OS in HCC patients, as indicated by a hazard ratio (HR) of 1.97, a 95% confidence interval (CI) of 1.51–2.6, and a P-value of less than 0.001, as determined through multivariate Cox analysis. It suggests that PCD-TME classifiers serve as independent prognostic markers for HCC patients. Furthermore, the PCD-TME classifiers demonstrated considerable predictive efficacy across various demographic and clinical factors, including gender, age, stage, and tumor grade, as depicted in Figures 5C–F. Within the ICGC-LIRI-JP dataset, the PCD-TME classifier was confirmed as a significant risk factor in hepatocellular carcinoma (HCC) (Figure 5G). Furthermore, we elucidated the correlation between the PCD-TME classifier and various clinical characters and validated the PCD-TME classifier in an independent cohort within this section.

This study investigated somatic alterations within the PCD-TME classifiers, with TCGA-LIHC identifying the 20 most frequent variant mutations (Figures 6A, B). Variations in TMB landscapes were observed between the two subgroups. TP53, TTN, CTNNB1, MUC16, and PCLO were among the top five mutations in the PCD high-TME low group, while CTNNB1, TTN, MUC16, ALB, and PCLO were among the top five in the PCD low-TME high group. Somatic mutations were more prevalent in the PCD low-TME high group compared to those shown in Figures 6A, B. Figure 6C demonstrates a notable distinction in TMB within the PCD-TME classifier. It has been established in previous research that CTNNB1 ranks among the proto-oncogenes exhibiting the highest mutation frequency in hepatocellular carcinomas (22%), with over half of hepatoblastomas presenting CTNNB1 mutations (Cleary et al., 2013; Harding et al., 2019; Nakagawa et al., 2019; Nakagawa and Shibata, 2013; Senni et al., 2019; Zucman-Rossi et al., 2015). The levels of CTNNB1 expression exhibited variability across the PCD-TME classifiers, as illustrated in Figure 6D. Notably, the subgroup characterized by PCD high-TME low exhibited a significantly elevated incidence of CTNNB1 mutation. Previous research has indicated that patients exhibiting low levels of CTNNB1 expression experienced longer survival durations in contrast to those with high CTNNB1 expression levels, which is consistent with our findings (Figure 6E). A comprehensive examination was conducted utilizing Kaplan-Meier curves in conjunction with CTNNB1 and PCD-TME subgroups. Notably, the classifier demonstrated the ability to discern patients with more favorable prognoses within the subset of individuals harboring CTNNB1 mutations (Figure 6F). These results suggest that the PCD-TME classifier exhibits heightened sensitivity in patient stratification and can effectively pinpoint improved prognostic outcomes in individuals with CTNNB1 mutations. Our hypothesis posited that the PCD-TME classifier could effectively forecast clinical responses in immunotherapy patients due to differences in immune statuses and tumor mutational burdens. To investigate this hypothesis, the TIDE algorithm was employed to predict responses to immunotherapy. Among patients with HCC who exhibited positive responses to immune checkpoint blockade (ICB) therapy, PCD scores were found to be significantly lower (Figure 6G). Additionally, the Proteomap tool was utilized to visually elucidate the potential mechanism through which the PCD-TME classifier predicts immunotherapy responses. As depicted in Figures 6H, I, Proteomap patterns exhibited striking similarities between individuals with PCD low-TME high levels, as well as those who responded favorably to immunotherapy. Furthermore, a notable resemblance was noted between the PCD high-TME low subgroup and the immunotherapy nonresponder (Supplementary Figure 2). In conclusion, these findings imply that a pretreatment PCD-TME signature could serve as a potential indicator of the patient’s tumor immune microenvironment, thereby aiding in the prediction of their therapeutic response.

Additionally, the validity of the PCD score was confirmed through single-cell transcriptomic analysis. Following quality control, a total of 19 cell clusters (Supplementary Figure 3) were identified from five samples. Subsequently, ten distinct cell types (Malignant, Macrophage, Endothelial, T/NK, Monocyte, Fibroblast, Plasma, DC, B cells, and Mast) were distinguished using marker genes (Figure 7A). The “AddModuleScore” function was utilized to compute PCD scores for individual cells, indicating that immune cells (Macrophages, Monocytes, and Mast Cells) displayed elevated PCD scores (Figures 7B, C). The examination of intercellular communication among immune cells was conducted to assess the potential impact of PCD status on their functionality. Macrophages were specifically chosen for further analysis based on the results depicted in Figure 7C. Figure 7D depicts the extent of interactions and the intensity of intercellular communication. Figure 7E shows that immune cells (macrophages) with a high expression of PCD are more likely to activate intercellular signaling with malignant cells. These findings highlight the pivotal role of PCD in modulating immune cell interactions, particularly macrophage-driven communication with malignant cells, suggesting PCD status may influence tumor-immune dynamics.

To further investigate potential PCD-related genes involved in regulating HCC development, the GEO dataset GSE10143 was utilized for cross-DEG analysis. Our findings reveal a significant upregulation of HTRA2 in HCC tissues compared to adjacent paracancerous, as illustrated in Figure 8A. Consistent results were achieved via WB, thereby validating the observed differential expression of HTRA2 within the overlapping DEGs (Figure 8B). The results of qRT-PCRT performed on L-O2 normal liver cells and JHH-7 HCC cells also show the similar results suggested by the crossed DEGs (Figure 8C). The HPA database was used to examine HTRA2 protein expression at the protein level. In the HPA database, there is a relatively high expression of HTRA2 in HCC when compared to normal (Figure 8D).