Introduction

Liver cancer ranks sixth in morbidity and third in mortality among all kinds of malignant tumors worldwide, with almost 906,000 new cases and 83,000 new deaths each year, making it a serious threat to human’s health.¹ Hepatocellular carcinoma (HCC) is the most common primary liver cancer and accounts for 80% of all cases. Its risk factors include hepatitis B virus (HBV) or hepatitis C virus (HCV) infection, alcohol ingestion, metabolic syndrome, and cirrhosis.^2–4 In terms of treatment, it is essential to develop an individualized and multidisciplinary approach based on the clinical stages and systemic conditions of specific patients with HCC. Liver resection and transplantation can be performed in patients with early-stage disease. For patients with unresectable advanced HCC, transarterial chemoembolization (TACE) and radiofrequency ablation (RFA) are performed.^5,6 Immunotherapy, systemic chemotherapy, and targeted therapy are widely used to treat advanced HCC. Several studies have reported that immunotherapy in combination with targeted therapy is superior to single-agent targeted therapy in terms of patient prognosis, which has become the preferred treatment for advanced HCC.^6–10 However, the median overall survival (OS) and progression-free survivals (PFS) are only 19.2 months and 6.8 months respectively for patients with HCC.¹¹ This may be attributed by the high heterogeneity of HCC. Therefore, it is necessary to differentiate patients with HCC more precisely at the molecular level.

Lysosomes are unique organelles that were first discovered by Prof. Christian de Duve in the 1950s.¹² Since their discovery, extensive research has revealed that lysosomes are not merely the cell’s recycling center but also play central roles in cellular degradation, metabolic regulation, and signal transduction.^13,14 They are critical for processes such as autophagy, immunoregulation, and nutrient sensing, and their dysfunction is implicated in both tumor and non-tumor diseases.^15,16 In cancers, lysosomes regulate tumor cell proliferation, invasion, and treatment resistance, making them a focal point for therapeutic exploration. For instance, Pechincha et al demonstrated that the lysosomal enzyme trafficking factor (LYSET) is essential for energy uptake in nutrient-deficient conditions, a common feature of the tumor microenvironment (TME).¹⁷ Notably, the growth of various tumors is significantly inhibited by LYSET knockout, highlighting the potential of targeting lysosomal pathways as a promising strategy for cancer treatment.¹⁸ These findings underscore the importance of investigating lysosome-related genes (LRGs) in HCC, as they may reveal novel molecular subtypes and provide valuable insights for personalized therapeutic approaches.

The TME is a dynamic and complex ecosystem that plays pivotal roles in cancer progression. Composed of immune cells (eg, macrophages, T lymphocytes, dendritic cells), cancer-associated fibroblasts (CAFs), vascular networks, the extracellular matrix (ECM), and a myriad of signaling molecules, the TME influences key cancer hallmarks such as apoptosis, evasion, proliferation, invasion, and metastasis.^19–22 Importantly, the TME has emerged as a promising therapeutic target, with numerous drugs now designed to modulate its components, including immune checkpoint inhibitors and anti-angiogenic agents.^23–27 Lysosomes are intricately linked to the TME, regulating critical biological processes such as the polarization of M2 macrophages, which promote tumor progression,²⁸ and the antigen-presenting capacity of dendritic cells, which shape anti-tumor immune responses.²⁹ Despite these advances, the specific mechanisms by which LRGs interact with the TME to drive HCC progression remain poorly understood. Elucidating these mechanisms is crucial, as it could uncover new therapeutic targets and pave the way for individualized treatment strategies for patients with HCC in the future.

In this study, bulk-RNA data for HCC and adjacent normal tissues were obtained from public databases. Differentially expressed LRGs were analyzed, and COX regression analysis was used to screen prognosis-related genes in HCC. Clustering analysis was used to identify two novel molecular subtypes based on the genes related to lysosomes in HCC. Furthermore, we analyzed pathway enrichment conditions, immune infiltration, immune checkpoint expression, and drug sensitivity in the two subtypes. In addition, we constructed a novel multiple-factors prognosis prediction model based on the screened LRGs. Finally, we hope that our results will provide new insights and guide clinical practice in the future.

Materials and Methods

Data Collection and Processing

Bulk-RNA data and clinical information of 374 patients with HCC were obtained from liver hepatocellular carcinoma (LIHC) cohort in The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). Data from the validation sets were downloaded from the International Cancer Genome Consortium (ICGC) database (https://docs.icgc-argo.org/docs/data-access/icgc-25k-data) and GSE database (GSE14520). 869 LRGs were acquired from Molecular Signature Database v7.5.1 (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/), which were shown in Supplementary Table 1. The entire analytical process of this study is presented in this flowchart (Figure 1).

Figure 1 Flow chart of this study.

Identification of Differentially Expressed LRGs

In this study, differentially expressed analysis was utilized to screen LRGs related to HCC with the “edgR” R package. To identify prognosis-associated LRGs, Cox proportional hazard regression analysis was performed using the prognosis data from the LIHC cohort of TCGA and the “survival” R package.

ConsensusClusterPlus Analysis of LRGs

The “ConsensusClusterPlus” R package was used to perform Consensus Cluster Plus analysis of clinical samples in TCGA-LIHC. The intersection of differentially expressed LRGs and prognosis-associated LRGs was selected as the key LRGs, which were mapped to HCC samples. By subsampling the proportion of features from these genes, an agglomerative hierarchical clustering algorithm was used to partition each sample into k clusters. The clustering process was repeated nine times. The clustering proportion of the key LRGs groups was defined by pairwise consensus values, which were calculated and integrated into the consensus matrix (CM), and the HCC samples were divided into k clusters. These clusters were defined as the consensus clusters (HCC subtypes). Moreover, the number of clusters and confidence in cluster sample qualifications were calculated based on quantitative and visual stability evidence. Principal component analysis (PCA) was also performed to minimize bias resulting from the categorical variable analysis.

Survival and Pathway Enrichment Analysis of Different Molecular Subtypes

Firstly, to compare the prognoses of the two subtypes, Kaplan-Meier (KM) survival analysis was performed using the “Kaplan-Meier survival” and “survminer” R packages. The Kaplan Meier survival curve was used to show the survival difference of the two subtypes. Then, with “Cluster profiler” R package, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) function and pathway enrichment analysis were performed to further explore downstream signaling pathways which key LRGs were mostly enriched in. Gene Set Variation Analysis (GSVA) was used to obtain differentially expressed gene sets in “h.all.v7.4. symbols.gmt” among the subtypes. Gene set enrichment analysis (GSEA) was also conducted on the selected genes expressing in both subtypes and the top 5 hallmark pathways were shown.

TME and Immune Infiltration Analysis

To analyze the differences of TME among different subtypes, Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) analysis was performed by utilizing the R package “estimate”. The infiltration degree of various immune cells in different subtypes was also calculated by Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm.^30–33 Single-Sample Gene Set Enrichment Analysis (ssGSEA) was used to evaluate and quantify the enrichment levels of immune-related sets in different subtypes.

Identification of Alternative Genes by Weighted Gene Correlation Network Analysis (WGCNA)

To further screen for alternative genes to construct a prognosis prediction model, WGCNA was performed. After determining the fittest soft-thresholding power for each subtype, the key LRGs were clustered into several modules. The degree of immune infiltration was considered a clinical phenotype for further co-expression analysis, and a module-trait relationship heatmap was created to find the module most associated with the degree of immune cell infiltration. The intersection of the results for each subtype was identified as an alternative gene used to construct the prognosis prediction model.

Independent Prognostic Analysis, Assessment of Clinical Relevance, and Construction of a Nomogram

To further assess the prognostic value of these alternative genes, a multivariate Cox regression analysis was performed. Subsequently, a prognostic model was constructed according to the formula (Risk Score = ), which could be utilized to predict the prognosis of patients with HCC. In this formula, “n” represented OS-associated LRGs, and “β” represented the coefficient of each OS-associated LRG. To further assess the clinical relevance of the prognostic model and OS-associated LRGs, gene expression heatmaps, survival status scatter plots, and risk-score distribution plots were constructed. Simultaneously, the AUC curves were used to test the accuracy of this prognostic mode using the R package “timeROC”. A nomogram was constructed to intuitively display the relationships between OS-associated LRGs and prognosis. The regression model was fitted and the nomogram was depicted utilizing the R package “rms”. The 1-year survival, 3-year survival, and 5-year survival rates could be illustrated clearly based on the expression of the selected key LRGs.

RT-qPCR Validation of Selected Key LRGs

Total RNA was isolated from ten pairs of HCC and normal tissues using TRIzol solution (Solarbio). The isolated RNA was reverse transcribed using a cDNA synthesis kit (Vazyme). RT-qPCR was performed using SYBR Green PCR Master Mix (Vazyme) and primers binding to the target genes, and GAPDH. The primer sequences were listed in Table 1. The cycling protocol was performed in accordance with the manufacturer’s instructions. The relative mRNA levels were calculated using the 2^−ΔΔCt method. These specimens were obtained from Liaoning Cancer Hospital and Institute, China. This study was approved by the Institutional Review Board of our institution.

Table 1 The Sequences of the Primers

Immune Checkpoints Analysis, Human Leukocyte Antigen (HLA) Analysis, and Tumor Immune Dysfunction and Exclusion (TIDE) Analysis

To study the immune mechanisms of each subtype, immune checkpoint and HLA analyses were performed. Using the TIDE analysis tool (http://tide.dfci.harvard.edu/), the immune escape phenomenon and overall response rate to immunotherapy in each subtype were also investigated.

Drug-Sensitive Analysis

The pRRophetic (version 0.5) algorithm in the R package was used to identify the most suitable targeted or chemotherapeutic drug. The half-maximal inhibitory concentration (IC₅₀) was calculated to assess the sensitivity of HCC cells to drugs using the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org).

Statistical Analysis

In this study, only a two-tailed P value < 0.05 was considered statistically different. R version 3.6.2 software and GraphPad Prism v. 8.01 were used for statistical analysis and visualizations. The Student’s t-test was performed to analyze the data obeying a normal distribution. The Mann–Whitney U-test was used to evaluate non-parametric data.

Results

Identification of Novel Molecular Subtypes of HCC Based on LRGs

The flowchart of the study was shown in Figure 1. A total of 109 genes were selected based on the interaction between the results of differential expression analysis and univariate Cox regression analysis (Figure 2A). These genes were used for the cluster analysis of 374 HCC samples. According to the consensus matrix plot, tSNE scatter plot, and PCA scatter plot (Figure 2B–D), HCC samples were divided into two subtypes, C1 and C2 (k=2). As shown in the heatmap, the expression of the selected LRGs in C1 subtype was significantly higher than that in C2 subtype (Figure 2E). In addition, the clinicopathological features of C1 and C2, including tumor size and pathological heterogeneity, were significantly different (Figure 2E). Therefore, C1 subtype could be regarded as a specific subtype with a high correlation with lysosome-related genes. The survival outcome of the C1 subtype was remarkably poorer than that of C2, according to the KM survival curve (Figure 2F). Subsequently, external validation was conducted using the ICGC and GSE databases and the same conclusions were obtained (Supplementary Figure 1). These results suggest that key LRGs are strongly linked to the prognosis of patients with HCC and deserve more in-depth analysis.

Figure 2 Identification of novel molecular subtypes of HCC based on LRGs. (A) The forest plot displayed the 109 key LRGs and their respective hazard ratios. (B) Cluster analysis by ConsensusClusterPlus. (C) Principal component analysis by tSNE. (D) Principal component analysis by PCA. (E) The heatmap showed HCC samples were divided into 2 subtypes, C1 subtype and C2 subtype based on the expression of the 109 key LRGs. (F) K-M survival curves showed different prognosis of the two HCC subtypes in TCGA.

Pathway Enrichment Analysis Between C1 and C2 Subtypes

Based on differential expression profiles, the two subtypes were significantly different at the molecular level. The differentially expressed genes were then used for pathway enrichment analysis of subtypes C1 and C2. We firstly performed GO analysis using TCGA, ICGC, and GSE databases and found that the immune response processes (immune response-regulating cell surface receptor signaling pathway involved in phagocytosis, inflammatory response, etc), positive regulation of cell communication, and positive regulation of signaling mainly differed between the two subtypes (Figure 3A–C). Moreover, GSEA-HALLMARK analysis was performed to avoid bias induced by the enrichment analysis. Using TCGA, ICGC, and GSE databases, HALLMARK analysis showed significant differences in the G2M_Checkpoint, E2F_Targets, Angiogenesis, KRAS signaling pathway, MTORC1 signaling pathway, and other important cancer-related pathways between the two subtypes (Figure 3D–F). Furthermore, using TCGA, ICGC, and GSE databases, we performed KEGG analysis and found that the two groups were mainly different in the cell cycle, ECM-receptor interaction, fructose and mannose metabolism, cellular senescence, proteoglycans in cancer, PPAR signaling pathway, fatty acid metabolism, p53 signaling pathway, biosynthesis of amino acids, ubiquinone, and other terpenoid-quinone biosynthesis pathways (Figure 3G–I). The results of the GSVA analysis also showed corresponding results (mTORC1_signaling, PI3K_AKT_mTOR_signaling, Wnt_beta-catenin_signaling, etc). (Supplementary Figure 2). These results suggest that the immune response, metabolic processes, and cancer-related signaling pathways play vital roles in shaping the two subtypes of HCC based on LRGs.

Figure 3 Pathway enrichment analysis between C1 and C2 subtypes. (A–C) Results of GO enrichment analysis based on TCGA, ICGC and GSE databases. (D–F) Results of GSEA-HALLMARK enrichment analysis based on TCGA, ICGC and GSE databases. (G–I) Results of KEGG pathway enrichment analysis based on TCGA, ICGC and GSE databases.

Immune Microenvironment Analysis of The Two Subtypes

To explore the differences in the immune microenvironment between subtypes C1 and C2, ESTIMATE, CIBERSORT algorithm, and ssGSEA analysis were further conducted. TCGA database was used as the test set, and the ICGC and GSE databases were used as validation sets. According to the ESTIMATE results, the immune Score of C1 subtype was significantly higher than that of C2 subtype (Figure 4A). Moreover, CIBERSORT and ssGSEA analyses showed that the proportions of immune cells between subtypes C1 and C2 were also different at the overall level (Figure 4B and C). Other cytokines, interferons and receptors, interleukins and receptors, chemokines and receptors involved in related intercellular communication pathways have also been investigated. The results suggested that the expression of cell communication molecules in C1 subtype was significantly higher than that in C2 subtype (Supplementary Figure 3A). The results were validated using ICGC and GSE databases (Supplementary Figure 3B and C). Communication molecules are strongly related to immune cells, suggesting that C1 subtype is closely related to immune cells, which also means that the degree of key LRGs expression in HCC patients may affect the efficacy of immunotherapy by regulating the level of immune infiltration.

Figure 4 Immune infiltration analysis of the two subtypes. (A) The box plots showed the stromal score, immune score and ESTIMATES score of both subtypes in TCGA, ICGC and GSE databases based on ESTIMATES. (B) The box plots showed percentage of several immune cells of both subtypes in TCGA, ICGC and GSE databases based on CIBERSORT algorithm. (C) The heatmaps showed the expression level of several immune cell gene sets of both subtypes in TCGA, ICGC and GSE databases based on ssGSEA. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001. Abbreviation: ns, non-significant.

Identification of Immune-Related LRGs by WCGNA

WGCNA was performed to further investigate the relationships between the key LRGs expression and immune cell infiltration. For the C1 subtype, the soft threshold was equal to 4, R² was equal to 0.87, and mean connectivity was equal to 6.16 (Figure 5A and B). Using the above parameters, a cluster dendrogram was constructed, showing key LRGs in the C1 subtype, which were clustered into three different gene modules marked by diverse colors (Figure 5C and D). The module-trait relationship heatmap illustrated that the turquoise module was strongly associated with different immune cells in the C1 subtype. The same method was used to analyze C2 subtype, with a soft threshold of 4. Then R² was equal to 0.90, and the mean connectivity was equal to 5.41 (Figure 5E and F). Key LRGs in C2 subtype can be clustered into four different gene modules, marked by diverse colors. The turquoise module in C2 subtype was relatively less related to immune cells (Figure 5G and H). Subsequently, the intersection of the turquoise module in subtypes C1 and C2 was performed, and alternative prognostic LRGs were identified for further analysis, including PRKCD, CPNE1, AGRN, P2RX4, PYGB, HRG, LAPTM4B, FTH1, SRC, DSN1, TUBB, LPCAT1, TTR, CKAP4, HSP90AB1, and GLA. The results also indicated that these identified genes were potentially critical players in the malignant development and prognostic prediction of the novel HCC subtypes C1 and C2, offering opportunities for deeper mechanistic insights, biomarkers identifying, and targeted therapeutic strategies development.

Figure 5 Identification of immune-related LRGs by WCGNA. (A–B) Analysis of network topology for soft powers in C1 subtype. (C) Dendrogram and genes module colors in C1 subtype. (D) The module-trait relationships displayed the correlations between different gene modules and immune cells in C1 subtype. (E–F) Analysis of network topology for soft powers in C2 subtype. (G) Dendrogram and genes module colors in C2 subtype. (H) The module-trait relationships displayed the correlations between different gene modules and immune cells in C2 subtype.

Construction of the Multiple Genes Prognosis Prediction Model

Multivariate Cox regression analysis was performed to reduce the dimensions of the input genes obtained from WGCNA analysis, PRKCD, LAPTM4B and LPCAT1 were selected as OS-associated LRGs and were used to construct a multiple gene prognosis prediction model. Using this model, all HCC samples were divided into high-risk and low-risk groups based on the median risk scores of all samples. Consistent with Figure 2, the high-risk group was similar to the C1 subtype and the low-risk group was similar to the C2 subtype (Figure 6A and B). The fishbone diagram also suggested that the risk score of C1 subtype was significantly higher than that of C2 subtype (Figure 6C). The heatmap showed that the expression of these three genes was remarkably higher in C1 subtype than that in subtype C2 (Figure 6D). The nomogram illustrates the survival outcomes of patients with HCC based on the expression of RKCD, LAPTM4B, and LPCAT1 in TCGA. With this help, the 1-year survival, 3-year survival, and 5-year survival rates could be intuitively predicted (Figure 6E). The AUC curve of ROC in TCGA was 0.730 at 1 year, 0.674 at 3 years and 0.660 at 5 years (Figure 6F); The AUC curve of ROC in ICGC was 0.778 at 1 year, 0.837 at 3 years and 0.729 at 5 years (Figure 6G); The AUC curve of ROC in GSE was 0.673 at 1 year, 0.672 at 3 years and 0.670 at 5 years (Figure 6H). The AUC values suggested that the model had good specificity and sensitivity.

Figure 6 Establishment of multiple factors prognosis prediction model. (A) Risk score distribution plots of C1 subtype and C2 subtype based on TCGA. (B) Survival status scatter plots of C1 subtype and C2 subtype based on TCGA. (C) The fishbone diagram showed the differences in risk score of both subtypes based on TCGA. (D) The heatmap showed the expression condition of OS-associated LRGs (PRKCD, LAPTM4B, and LPCAT1) in both subtypes. (E) Nomogram developed based on the expression of PRKCD, LAPTM4B, and LPCAT1. (F) ROC curves for OS-related LRGs based on TCGA database. (G) ROC curves for OS-related LRGs based on ICGC database. (H) ROC curves for OS-related LRGs based on GSE database. ***P<0.001.

The Expression of Selected Three Genes Were Significantly Increased in HCC Tissues Compared with Normal Tissues

We further validated the expression of three selected genes, PRKCD, LAPTM4B, and LPCAT1 in ten pairs of fresh HCC and normal tissues using RT-qPCR analysis. We found that the expression levels of PRKCD, LAPTM4B, and LPCAT1 were significantly higher in HCC tissues than those in normal tissues (P<0.05, Figure 7A–C). These results indicated that they may serve as potential therapeutic targets for patients with HCC.

Figure 7 Validation of the high expression of PRKCD (A), LAPTM4B (B), and LPCAT1 (C) in HCC and normal samples via RT-qPCR analysis. *P<0.05.

Immune Checkpoints and TIDE Analysis

Furthermore, we explored the landscape of immune checkpoints for both subtypes. Compared to the C2 subtype, all types of immune checkpoints and HLA were highly and ubiquitously expressed in the C1 subtype using the data (Figure 8A and B), which was also validated in the ICGC and GSE databases (Figure 8C–F). Meanwhile, the bar chart and circular chart suggested that the immunologic response rate in C1 subtype was significantly lower than that in C2 subtype based on data from the ICGC, and GSE databases, expect the data in TCGA databases (Supplementary Figure 4). The violin chart also showed that the TIDE score of the C1 subtype was significantly higher than that of the C2 subtype, suggesting that the possibility of the immune escape phenomenon in the C1 subtype may be higher than that of the C2 subtype based on the data from the above databases (Figure 8G–I).

Figure 8 Immune checkpoints, HLA and TIDE analysis of the two subtypes. (A) The results of immune checkpoints analysis in both subtypes based on TCGA. (B) The results of HLA analysis in both subtypes based on TCGA. (C) The results of immune checkpoints analysis in both subtypes based on ICGC. (D) The results of HLA analysis in both subtypes based on ICGC. (E) The results of immune checkpoints analysis in both subtypes based on GSE. (F) The results of HLA analysis in both subtypes based on GSE. (G–I) The results of TIDE analysis in both subtypes based on TCGA, ICGC and GSE databases. *P<0.05, **P<0.01, ***P<0.001. Abbreviation: ns, non-significant.

Drug-Sensitive Analysis in Multiple Databases

Based on data from the GDSC website, drug sensitivity responsiveness was also assessed between subtypes C1 and C2. We found that the responsiveness of candidate drugs between subtypes C1 and C2 was relatively different and that C1 subtype had higher drug sensitivity responsiveness. We found that the IC₅₀ values of the Akt inhibitor A.443654, PLK1 inhibitor BI.2536, Gemcitabine, SGK1 inhibitor GSK.650394, Pyrimethamine, Eg5 inhibitor S.Trityl.L.cysteine in C1 subtype were much lower in the C1 subtype. These drugs might provide more efficient treatments for patients with HCC (Figure 9A–C).

Figure 9 Drug sensitivity analysis. (A) The results of drug sensitivity analysis in both subtypes based on TCGA. (B) The results of drug sensitivity analysis in both subtypes based on ICGC. (C) The results of drug sensitivity analysis in both subtypes based on GSE.

Discussion

HCC is the most common primary liver cancer,^34–36 and is treated with multidisciplinary management, including hepatectomy, liver transplantation, radiofrequency ablation, chemotherapy, immunotherapy, and targeted therapy at present.^37–40 However, tumor recurrence, metastasis, and chemotherapy resistance greatly affect the survival outcomes of patients with HCC.^41,42 Researchers have constructed several molecular subtypes to guide HCC therapies and improve the prognosis of patients with HCC. The molecular subtypes include metabolism, immunity, polyposis, cuproptosis, ferroptosis, glycolysis, and so on.^43–50 Although various molecular subtypes have been identified, the survival status of patients is still not very satisfied. Thus, the development of novel and accurate molecular subtypes to improve HCC survival is urgently required.

Lysosomes play an important role in various diseases including various kinds of cancers. Zeng et al reported that lysosomal degradation was targetable in osteosarcoma lung metastasis.⁵¹ Alejandro et al found that the lysosomal pathway contributed to cisplatin-induced cancer cell death.⁵² Therefore, it may be useful to predict and improve the survival outcomes of patients with HCC by constructing novel molecular subtypes based on the LRGs. In this study, 869 LRGs were downloaded from MSigDB, and differential expression analysis was performed based on the Bulk-RNA sequence data of HCC and normal tissues from TCGA database. After differential expression and survival analyses, 109 screened genes were used for cluster analysis. Using consensus clustering analysis, 374 HCC samples were clustered into two subtypes based on the 109 prognostic LRGs. Remarkably, the survival time in the C1 subtype was significantly shorter than that in the C2 subtype, indicating that the LRG-related molecular type could significantly distinguish the survival outcomes of patients with HCC. The prognostic differences between subtypes C1 and C2 also suggest that LRGs play a vital role in HCC. Thus, we further conducted GO, GSEA, KEGG, and GSVA enrichment analyses to explore the differences between the downstream signaling pathways of both subtypes. The differentially modulated genes were significantly linked to cancer-related pathways, including hallmark_G2M_checkpoint, hallmark_E2F_targets, p53 signaling, mTORC1 signaling, PI3K_AKT_mTOR_signaling, Wnt_beta_catenin_signaling, Notch signaling, TNFA_signaling_via_NFκB, and the inflammatory_response. In addition, glycolysis, fructose and mannose metabolism, cholesterol homeostasis, and fatty acid metabolism related to solid cancers were also altered in the enrichment analysis. Moreover, previous studies have reported that hallmark_G2M_checkpoint, hallmark_E2F_targets, mTORC_signaling pathway, and PI3K_Akt_mTOR signaling pathway are strongly related to tumor proliferation and metastasis.^53–56 Glycolysis, cholesterol homeostasis, and fatty acid metabolism were also closely associated with the metastasis and invasion of solid tumors.^57–59 In addition, the pathways related to anti-tumor immune responses, including leukocyte migration, phagocytosis, inflammatory response,⁶⁰ were also significantly enriched between these two subtypes, suggesting that the two subtypes might have differences in the infiltration degree of several immune cells. The following analyses further support this hypothesis.

Tumor occurrence, invasion, metastasis, and progression are significantly affected by tumor microenvironment components such as various immune cells, stromal cells, and cytokines.^61,62 Immune cells play an important role in this process.^63,64 Studies have reported that Tregs act as anti-inflammatory cells, inhibit the immune response, and promote the malignant progression of tumors by secreting anti-inflammatory cytokines.⁶⁵ Neutrophils also play a vital role in accelerating tumor occurrence and metastasis by secreting anti-inflammatory cytokines to regulate TME.⁶⁶ In addition, cytotoxic immune cells, such as CD8⁺ and γδ T cells, secrete pro-inflammatory cytokines to exert anti-tumor effects.^67,68 Therefore, immune infiltration analysis was performed in HCC using ssGSEA, CIBERSORT, ESTIMATE, and WGCNA methods. M0 macrophages, Tregs, dendritic cells, and neutrophils were present in a higher proportion in the C1 subtype, whereas CD8⁺ T cells and γδ T cells had a relatively lower percentage in the C1 subtype. These results indicated that compared with C2 subtype, C1 presented an immunoinhibitory phenotype, which was consistent with the worse survival outcomes of patients in C1 subtype.

Using WGCNA, we further screened for immune-related LRGs at the intersection of immune-related modules in both subtypes. Multivariate Cox regression analysis was performed based on the screened LRGs, and three OS-associated LRGs were identified. OS-associated LRGs include lysosomal protein transmembrane 4 beta (LAPTM4B), protein kinase C delta (PRKCD), and lysophosphatidylcholine acyltransferase 1 (LPCAT1). The multifactor prognosis model has good prediction performance and an excellent ability to distinguish the prognosis of patients with HCC. Prognostic scores can also be visually acquired using nomogram. LAPTM4B is regarded as a diagnostic biomarker and potential therapeutic target for HCC. Wang et al found that LAPTM4B significantly promoted the proliferation and autophagy of HCC.⁶⁹ Wang et al found that LAPTM4B contributed to malignant proliferation of HCC stem cells and migration of myeloid-derived suppressor cells (MDSCs).⁷⁰ Meanwhile, under the positive action of the transcription factor AP4, LAPTM4B remarkably promoted the invasion and metastasis of HCC and reduced sensitivity to chemotherapy.⁷¹ Moreover, using clinical samples, Zhai et al identified that allele*2 of LAPTM4B may be significantly related to the genetic susceptibility to HCC and is considered a prognostic risk factor for patients with HCC.⁷² Li et al found that PRKCD could serve as a key messenger in the interaction between HCC cells and platelets and play a crucial role in the process of platelet-induced tumor progression.⁷³ Xu et al also found that high expression of PRKCD predicted poor survival outcomes in patients with HCC and was involved in tumor immune escape.⁷⁴ LPCAT1 also contributes to malignant development of HCC. Sun et al found that silencing of LPCAT1 remarkably inhibited the proliferation and metastasis of HCC by targeting S100A11 and Snail.⁷⁵ Zhang H and He RQ et al also found upregulated expression of LPCAT1 in HCC samples compared to non-tumor samples, which may serve as an independent prognostic predictor for the survival ofpatients with HCC. In their research, knocking down of LPCAT1 remarkably inhibited the proliferation, migration, and metastasis of HCC cells.^76,77 In our study, we validated the high expression of PRKCD, LAPTM4B, and LPCAT1 in 10 pairs of fresh HCC and normal samples using RT-qPCR analysis. All results indicate that our nomogram may have a relatively high predictive ability for the prognosis of patients with HCC.

In addition to the level of immune cell infiltration, the influence of immune checkpoints on the prognosis of patients with HCC cannot be neglected. Immune checkpoints are a crucial part of immune reactions and can inhibit certain anti-tumor immune responses, which are indispensable for regulating and maintaining immune homeostasis.^78–80 For example, immune checkpoint B and T lymphocyte attenuators (BTLA) can cause immunosuppression by inhibiting the activation and proliferation of B and T cells.⁸¹ A previous study also reported that high expression of BTLA was positively related to the expression of PD-L1 and indicated poorer prognosis in patients with non-small cell lung cancer.⁸² In addition, classical Cytotoxic T lymphocyte associated protein-4 (CTLA-4) plays a vital role in inhibiting T cell activation and serves as a biomarker for poor prognosis in several solid cancers including HCC,^83–85 nasopharyngeal carcinoma,⁸⁶ breast cancer,⁸⁷ and osteosarcoma.⁸⁸ In our study, we found that the expression of various immune checkpoints genes in C1 subtype was significantly higher than that in C2 subtype, such as BTLA, lymphocyte activating 3(LAG-3), CTLA-4 and so on. Meanwhile, based on the TIDE analysis, the score of C1 subtype was much higher than that of C2 subtype, which meant that immune evasion was more likely to happen in C1 subtype. Moreover, the higher proportion of non-responders in the C1 subtype indicated that the efficiency of immunotherapy in patients with HCC was poorer. Thus, the above results might explain why patients with HCC in the C1 subtype had a worse survival outcome. In addition, for patients who cannot undergo curative resection, systemic chemotherapy and targeted therapy are alternative methods.^89–92 To explore the drug sensitivity of both subtypes, drug-sensitivity analysis was performed. We found that A.443654, BI.2536, Gemcitabine, GSK.650394, pyrimethamine, and S. Trityl. L. cysteine were more sensitive in patients with C1 subtype than with C2 subtype. These drugs may serve as drug candidates for treating patients with C1 subtypes and deserve further research.

This study also has some limitations. Firstly, data from public databases were limited, and the clinical information of some samples might miss, resulting in biases and errors. Secondly, our conclusion was mainly obtained from bioinformatics, and more validations are needed via in vivo and in vitro assays. Thirdly, this study was conducted based on retrospective data, and it is necessary to design a prospective clinical trial with a large sample size to prove our hypothesis.

Conclusions

In conclusion, our study revealed the heterogeneity of HCC. We not only defined two new molecular subtypes of HCC based on differentially expressed LRGs, but also studied the correlations between molecular subtypes and prognosis by bioinformatic analysis of bulk sequencing data. We further explored the differences in the signaling pathways, immune cell infiltration, and sensitivity to drug therapy between the two subtypes. Additionally, we constructed a novel multiple-factor prognosis prediction model based on the screened LRGs. Finally, we hope the LRG-related subtypes of HCC will provide new insights and guidelines for future clinical practice.