System analysis based on the lysosome-related genes identifies HPS4 as a novel therapy target for liver hepatocellular carcinoma

Background: Liver cancer is a leading cause of cancer-related deaths worldwide. Lysosomal dysfunction is implicated in cancer progression; however, prognostic prediction models based on lysosome-related genes (LRGs) are lacking in liver cancer. This study aimed to establish an LRG-based model to improve prognosis prediction and explore potential therapeutic targets in liver cancer.

Methods: Expression profiles of 61 LRGs were analyzed in The Cancer Genome Atlas liver cancer cohorts. There were 14 LRGs identified, and their association with clinical outcomes was evaluated. Unsupervised clustering, Cox regression, and functional assays were performed.

Results: Patients were classified into high-risk and low-risk subgroups based on the 14 LRGs. The high-risk group had significantly worse overall survival. Aberrant immune infiltration and checkpoint expression were observed in the high-risk group. Furthermore, HPS4 was identified as an independent prognostic indicator. Knockdown of HPS4 suppressed liver cancer cell proliferation and induced apoptosis.

Conclusion: This study developed an LRG-based prognostic model to improve risk stratification in liver cancer. The potential value of HPS4 as a therapeutic target and biomarker was demonstrated. Regulation of HPS4 may offer novel strategies for precision treatment in liver cancer patients.

1 Introduction

Globally, LIHC ranks sixth among the most common cancers (1, 2). There are projected to be more than one million cases of liver cancer by 2030 (3). The number of patients who die from LIHC each year exceeds half a million (4). In the early stages, LIHC has no symptoms but rapidly progresses (5). At present, LIHC is mainly treated with liver transplantation, hepatic resection, and medication (6). LIHC can be surgically excised in only 10%–20% of patients, but recurrences are common (7). The current drug therapy has limited efficacy for LIHC because it is chemoresistant (8, 9). Despite the application of new treatment strategies for LIHC, efficacy remains unsatisfactory (10). Thus, identification of specific prognostic markers is essential for guiding LIHC therapy and improving OS of patients.

In all types of cells, lysosomes are membrane-bound phospholipid bilayers that catabolize protein degradation and recycle it through phagocytosis, endocytosis, and autophagy (11–13). It has been reported that lysosome functions are drastically altered during cancer progression, including alterations in lysosome volume, localization, and composition within the cell (14, 15). In addition, lysosomes have been found to be routed to the periphery of the cell under multiple stimuli that can result in a variety of pathological conditions (11). In malignant transformations, lysosomes are juxtaposed to the plasma membrane, which contributes to cell invasion and migration (16, 17). A number of cancers show significant overexpression of lysosomal hydrolases, which increases invasion of tumor. TMEM106B has been reported to increase lysosomal hydrolase synthesis in lung cancer cells, which is packaged into lysosomes and enlarge lysosomes (11). A hyperinvasive microenvironment is created as lysosomes secrete their protease cargo into the extracellular matrix under conditions of calcium flux (18–20).

We identified 14 genes out of 61 lysosome-related genes (LRGs) as significant predictors in LIHC. The 14 genes are ANKRD27, AP3M1, BCL10, CD63, CTSC, GLA, HPS1, HPS4, NPC2, PPT1, RAB7A, TPP1, VPS45, and RAMP3. The 14 lysosome-related genes were screened to classify molecular subgroups by NMF consensus clustering. The patients of LIHC from TCGA data were divided into C1 and C2 clusters. In addition, we explored that the most common CNVs are heterozygous amplifications and deletions of genes based on the LRGs. Further analysis of the pathway prediction and methylation was performed. An OS prognostic signature was constructed using LASSO regression analysis based on the LRGs. As an external verification database, we used the ICGC-LIHC cohorts for training the prognosis model.

We further explored the association between LRGs, immune checkpoints, and immune cells. Moreover, the TIDE scores of group G1 were significantly higher than those of group G2. The chemotherapy response in the two subtypes was also explored by comparing the IC50s of doxorubicin and cisplatin. We developed the nomogram by the univariate Cox regression model. Our result demonstrated that HPS4 and a worse prognosis are strongly associated with LIHC. We further explored the relation between HPS4 and immune cells based on eight single-cell datasets. We found that HPS4 could reduce proliferation and apoptosis significantly in liver cancer cells; it may be a promising therapeutic target and biomarker for LIHC, and regulating its activity may be an effective treatment strategy.

2 Materials and methods

2.1 Differential expression analysis of LRGs

We obtained 61 LRGs from the MSigDB (Supplementary Table 1). Expression profiles and clinical information for LIHC are provided via TCGA dataset (Supplementary Table 2). The raw read counts of genes were converted to transcripts per million (TPM) using STAR. Further analysis of LRG expression in LIHC was carried out using version R 3.6.3 and ggplot2. Our statistical significance criterion of log2(fold change) > 2 and adjusted p-value <0.05 indicated a significant difference.

2.2 Nomogram construction

Our goal was to identify the right terms for the nomogram by using univariate Cox regression analysis. For each variable, we calculated p values, HRs, and 95% confidence intervals using the “forestplot” R package. In order to predict 1-, 3-, and 5-year overall recurrence, a nomogram was developed using univariate Cox proportional hazards analysis. An individual’s recurrence risk can be calculated using nomograms by the “rms” R package.

2.3 The association between LRGs and gene mutation, CNV

Data about somatic mutations were provided by the UCSC Xena server for the GDC TCGA-LIHC project. According to the mutation order, the result was generated using the R package “maftools”. Data from TCGA database were downloaded and processed using GISTIC2.0, which can identify regions significantly altered including amplification and deletion.

2.4 The analysis of correlation between LRGs and immune infiltration

A TCGA dataset was downloaded to obtain LIHC clinical information and RNA sequencing profiles (level 3). LRG gene expression was correlated with immune scores using R’s ggstatsplot package, and multigene correlations were analyzed using R’s pheatmap package. To describe correlations between quantitative variables that do not follow a normal distribution, Spearman’s correlation analysis was used. R was used to implement all the analysis methods, and a statistical significance level below 0.05 was considered statistically significant.

2.5 Construct prognostic signature of LRGs

A TCGA dataset was downloaded to obtain LIHC clinical information and RNA sequencing profiles (level 3). In addition to converting count data to TPM, log2 (TPM+1) is normalized; it is important to keep samples together that contain information about the patient’s clinical status. A log-rank test was used to compare survival rates between the two groups. TimeROC (v0.4) was used to compare predictive accuracy between LRGs and risk scores. The features were selected using LASSO regression, cross-validated 10 times, and analyzed with the R package glmnet. With the survival package and Cox regression analysis, a prognostic model was constructed. We generated Kaplan–Meier curves with log-rank p-values and 95% confidence intervals (CIs) using log-rank tests. R was used to implement all the analysis methods, and a statistical significance level below 0.05 was considered statistically significant.

2.6 Identification of potential subtypes

80% of the samples were analyzed 100 times by ConsensusClusterPlus R package (v1.54.0). Using the R software package pheatmap (v1.0.12), we generated cluster heatmaps. Heatmaps showing gene expression were retained for genes with SD over 0.1. When the input gene number exceeded 1,000 after sorting the SD, the top 25% were extracted. All analysis methods and R packages were implemented in R version 4.0.3. 80% of TCGA LIHC samples (n = 371) were used as a discovery cohort to identify molecular subgroups. The remaining 20% of samples (n = 93) were held out as an internal validation set, to evaluate the reproducibility of the subgroups. Consensus clustering was run on the 80% discovery cohort, identifying two stable subgroups. These two subgroups were confirmed by testing their ability to classify the held-out 20% of samples.

2.7 The correlation between subtypes and immune

Evaluation of immune scores was conducted using immunedeconv. We extracted the expression of immune checkpoints including HAVCR2, PDCD1, CTLA4, LAG3 SIGLEC15, TIGIT, CD274, and PDCD1LG2. We predicted potential ICB responses using the TIDE algorithm and calculated mRNAsi based on the OCLR algorithm (21). Based on the mRNA expression signature, 11,774 genes are present in the gene expression profile. In order to map the dryness index to the range, the minimum value was subtracted and divided by the maximum value. R was used to implement all the analysis methods, and a statistical significance level below 0.05 was considered statistically significant.

2.8 The association of LRGs and drug sensitivity

With the help of the Genomics of Drug Sensitivity in Cancer (GDSC), we were able to predict chemotherapeutic response based on each sample. A prediction process was implemented using R package “pRRophetic”. The half-maximal inhibitory concentration (IC50) of samples was calculated using ridge regression. Using the batch effect of battle and tissue type, the duplicate gene expression values were summed up as a mean. R was used to implement all the analysis methods, and a statistical significance level below 0.05 was considered statistically significant.

2.9 Cell culture and transfection

From Procell (Wuhan, China), SK-hep-1 and Huh7 cells were obtained without mycoplasma infection. At 37°C with 5% CO₂, the SK-hep-1 and Huh7 cell lines were cultured in MEM (Gibco, Grand Island, NY) and DMEM (Gibco), respectively. The medium includes 10% fetal bovine serum (FBS; Gibco), 100 U/mL penicillin, and 100 mg/mL streptomycin (Invitrogen, Waltham, MA). We transfected si-HPS4 and si-NC using Lipofectamine 2000 (Invitrogen). The cells were harvested 24 h after transfection.

2.10 Cell proliferation and apoptosis

Cells transfected with si-HPS4 were seeded in cell culture plates, 5,000 cells per well. A colony formation assay was performed by inoculating 500 cells into six-well plates and growing them at 37°C for 15 days. Following fixing the colonies with methanol and staining them with hematoxylin, TRIzol (Invitrogen) was used to extract RNA from liver cancer cells. All of the primer sequences and small interfering RNA sequences are listed in Supplementary Table 2. A 48- h transfection of si-HPS4 or si-NC was followed by harvesting of SK-hep-1 and Huh7 cells. The cells were then stained with FITC Annexin V Apoptosis Detection Kit (Meilunbio, Dalian, China) and detected using flow cytometry (Beckman Coulter, Brea, CA). We analyzed the results using the FCAP Array software. A triplicate of each trial was conducted.

2.11 Single-cell RNA sequencing-based analysis

We further explored the correlation between HPS4 and immune cells based on eight single-cell datasets by the TISCH2 resource (22). We next analyzed the expression of HPS4 in the hepatocyte population and hepatocellular carcinoma cells by single-cell and HCL resources (23). We further explored the gene expression levels in all cancer single-cell samples by the CancerSCEM resource (24).

3 Results

3.1 Identification of prognostic molecular subgroups

There were 14 genes out of 61 lysosome-related genes (LRGs) identified as significant predictors in LIHC. Among them, 13 are risk factors, namely, ANKRD27, AP3M1, BCL10, CD63, CTSC, GLA, HPS1, HPS4, NPC2, PPT1, RAB7A, TPP1, and VPS45, whereas RAMP3 is a favorable factor (Figure 1A). According to their cumulative distribution function and function delta area, k = 2 appeared to be the best clustering value for the 14 genes (Figure 1B). Initially, the 14 lysosome-related genes were screened to classify molecular subgroups by NMF consensus clustering. The patients of LIHC from TCGA data were divided into C1 and C2 clusters according to the consensus map (Figure 1C). Heat maps of differentially expressed genes show red for high expression and blue for low expression (Figure 1D). An additional confirmation was obtained by performing a principal component analysis (PCA) (Figure 1E). The log-rank test was used to examine the survival of the different groups by utilizing Kaplan–Meier curves (p < 0.0001) (Figure 1F). We found that cluster 2 had a worse clinical prognosis compared with cluster 1. Our result demonstrated that LRGs were associated with drug sensitivity based on the CTRP and GDSC resources (Supplementary Figure 1). The differential genes between cluster G1 and G2 were further explored (Supplementary Figure 3).

FIGURE 1

Figure 1 Clustering of molecular subgroup based on the lysosome-related genes (LRGs) in liver hepatocellular carcinoma (LIHC). (A) There were 14 genes out of 61 LRGs identified as significant predictors in LIHC. (B) K = 2 appeared to be the best clustering value for 14 genes by cumulative distribution function and function delta area. (C) There were 14 lysosome-related genes screened to classify molecular subgroups by NMF consensus clustering. (D) Heat maps of differentially expressed genes between the C1 and C2 clusters. (E) Principal component analysis (PCA). (F) Cluster 2 had a worse clinical prognosis compared with cluster 1. Statistical analyses were conducted by one-way ANOVA, principal component analysis, and log-rank test. p < 0.05 was considered statistically significant.

3.2 Analysis of LRG copy number variations and immune correlations

A comparison of LRG expression between normal and tumor tissues was performed in LIHC. The LRG expression in tumors was obviously higher than in normal tissues, whereas RAMP3 expression was significantly lower in tumors (Figures 2A, C). A strong correlation was observed according to LRG expression, such as RAB7A and AP3M1 which are most correlated (r = 0.72) (Figure 2B). There was a significant positive correlation between CNV and mRNA expression in most LRGs (Figure 2D). The copy number variation (CNV) is prevalent in human cancers, contributing to tumor development. We found that the most common CNVs are heterozygous amplifications and deletions of genes, whereas homozygous amplifications and deletions are much less common (Figures 2E, G, H). The impact of tumor CNV levels on immune evasion is particularly intriguing. Analyzing the differences in immune infiltration between gene-set CNV groups, we estimate the relationship between immune and gene- set CNV levels. In the amplification group, NK cells, CD4 T cells, Tfh cells, Th2 cells, CD4 naïve cells, and macrophage, their infiltration scores all showed a negative correlation with CNV, whereas that of Treg cells was positively correlated with CNV in the group of amplification. The macrophage, Tfh cell, and CD4 naïve cell infiltration scores were negatively correlated with CNV in the group of deletion; however, the B- cell score showed the opposite trend (Figure 2F).

FIGURE 2

Figure 2 The copy number variation (CNV) of LRGs was correlated with immune cells. (A, C) A comparison of LRG expression between normal and tumor tissues was performed in LIHC. (B) A strong correlation was observed according to the expression of LRGs. (D) There was a significant positive correlation between CNV and most LRGs. (E, G, H) The most common CNVs are heterozygous amplifications and deletions of genes, whereas homozygous amplifications and deletions are much less common based on the LRGs’ (F) The relationship between immune and gene set CNV level was estimated. Statistical analyses were performed using t-test, Pearson correlation, and Fisher’s exact test. p < 0.05 was considered statistically significant. **p < 0.01; ***p < 0.001.

3.3 Pathway enrichment analysis and DNA methylation of LRGs

Further analysis of the CRGs’ GSVA (gene set variation analysis) scores was performed. The GSVA is an unsupervised method for measuring variations in gene- set activity (represented by GSVA scores) across cancer populations. By GSVA, the activity of differential gene sets between tumor and normal samples is analyzed. We found that the GSVA score was significantly increased in LIHC (Figure 3A). There was a worse prognosis in LIHC for patients with a high GSVA score (Figure 3B). In addition to the GSVA score, we studied the correlation between tumor-related pathways in the LIHC and the GSVA score. Our results demonstrated that GSVA score was positively associated with apoptosis and EMT pathway, and GSVA score was negatively associated with hormone ER and RTK pathway (Figure 3C). Additionally, LRG expression was associated with immune cell infiltrates using Spearman’s correlation analysis. We found that LRGs were negatively associated with Th17 cell, neutrophil, and monocyte, whereas LRGs were positively associated with cytotoxic cell, NK cell, Tfh cell, Treg cell, macrophage, CD8 T cell, Th1 cell, NKT cell, and DC cell (Figure 3D). In addition, we found that PPT1, HPS1, and RAMP3 were hypermethylated and CTSC, CD63, and GLA have a hypomethylation level (Figure 3E). Methylation and mRNA expressions were also investigated. We explored that most of LRGs had a negative association with methylation including ANKRD27, RAB7A, AP3M1, BCL10, NPC2, CTSC, HPS4, VPS45, HPS1, PPT1, and CD63 whereas GLA had a positive association with methylation (Figure 3F). In the LIHC cohorts, hypomethylation of VPS45 had an association with a worse prognosis whereas hypermethylation of GLA showed the opposite situation (Figure 3G). Our results demonstrated that LRGs had potential impact on pathways including the TSC/mTOR, RTK, RAS/MAPK, PI3K/AKT, hormone ER, hormone AR, EMT, DNA damage response, cell cycle, and apoptosis pathways. Most LRGs could activate the apoptotic cell cycle and the EMT pathway, inhibiting the RTK and TSC/mTOR pathways (Figures 3I, H). Immunohistochemistry staining images from the HPA were analyzed to further explore the LRG protein expression in LIHC. Our results demonstrated that the expression of most LRG proteins was higher in tumor tissues (Figure 4).

FIGURE 3

Figure 3 The analysis of pathway prediction and methylation. (A) The activity of differential gene sets between tumor and normal samples is analyzed based on the GSVA (gene set variation analysis) scores. (B) There was a worse prognosis in LIHC for patients with a high GSVA score. (C) The correlation between tumor-related pathways in LIHC and GSVA score was analyzed. (D) LRGs were associated with immune cell infiltrates using Spearman’s correlation analysis. (E) The hypomethylation level was explored based on the LRGs. (F) The association between methylation and mRNA expression was investigated. (G) The hypomethylation level of LRGs was associated with clinical prognosis. (H, I) The LRGs had a potential impact on pathway. Statistical analyses included t-test, Pearson correlation, log-rank test, and enrichment analyses. p < 0.05 was considered statistically significant. *p < 0.05; #FDR ≤ 0.05.

FIGURE 4

Figure 4 Immunohistochemistry staining images from the HPA resource were analyzed to further explore the LRGs protein expression in LIHC.

3.4 Development of an LRG-based prognostic signature model

As part of this research, we evaluated the relationship between LRGs and the prognosis of LIHC patients. Based on univariate Cox regression, LRGs were selected for LASSO regression (Figures 5A, C). An OS prognostic signature was constructed using LASSO regression analysis for 10 genes of LRGs: risk score = (0.0817)*BCL10+(0.0043)*CD63+(0.0601)*CTSC+(0.0502)*GLA+(−0.2469)*HPS4+(0.1112)*PPT1+(0.3027)*RAB7A+(-0.2508)*RAMP3+(0.1498)*TPP1+(0.1685)*VPS45 (lambda.min = 0.0173). There was a significant association between LRGs’ signature risk score and poor overall survival in the study (HR = 2.312 (1.619–3.301), log-rank p = 3.99e−06). In terms of AUCs, the 1-, 3-, and 5-year ROC curves were accurate by 0.774, 0.737, and 0.723, respectively. Survival of LIHC patients was significantly associated with LRG-related risk signatures at the individual level (Figure 5E). As an external verification database, we used the ICGC-LIHC cohorts for training the prognosis model. Based on univariate Cox regression, LRGs were selected for LASSO regression in the ICGC-LIHC cohorts (Figures 5B, D). There was a significant association between LRGs’ signature risk score and poor overall survival in the ICGC-LIHC cohorts (HR = 3.696 (1.821–7.501), log-rank p = 0.000296, Figure 5F). The association between subtypes and ferroptosis and m6A was further performed (Supplementary Figure 2). We also found that the relative expression levels of macrophages were lower in the high-risk group (Supplementary Figure 4).

FIGURE 5

Figure 5 Prognostic signature construction for LRGs. (A, C) Based on univariate Cox regression, LRGs were selected for LASSO regression. (B, D) As an external verification database, we selected the ICGC-LIHC cohorts for training the prognosis model. (E) Survival of LIHC patients was significantly associated with LRG-related risk signatures at the individual level in the TCGA-LIHC cohorts. (F) There was a significant association between LRGs’ signature risk score and poor overall survival in the ICGC-LIHC cohorts. Statistical analyses utilized LASSO regression, log-rank test, and time-dependent ROC analysis. p < 0.05 was considered statistically significant.

3.5 Evaluation of subtype-specific immune landscape and drug response

A comparison of LRG expression in the two subtypes was also carried out. We found that most LRGs increased significantly in group G1 whereas gene RAMP3 showed the opposite trend (Figure 6A). In tumors with higher mRNAsi, the tumor has dedifferentiated more and its cancer stem cells are more active. With a better clinical prognosis, the G2 group had lower mRNAsi scores whereas the G1 group had higher mRNAsi scores (Figure 6B). Our results demonstrated that ANKRD27, AP3M1, BCL10, CTSC, HPS1, HPS4, PPT1, RAB7A, TPP1, and RAMP3 had a negative relation to immune infiltration of CD4 T cells using the EPIC algorithm. Moreover, ANKRD27, AP3M1, BCL10, CD63, CTSC, GLA, HPS1, HPS4, NPC2, PPT1, RAB7A, TPP1, and VPS45 had a positive correlation with macrophage (Figure 6C). Furthermore, the TIDE algorithm showed low efficacy of immune checkpoint blockade therapy (ICB) based on high TIDE scores. As a result, the TIDE scores of group G1 were significantly higher than those of group G2 (Figure 6D). We found that subtypes had a significant association with immune cells including B cell, CD4 T cell, CD8 T cell, endothelial cell, macrophage, and NK cell. Macrophages in the G1 group had lower immune scores than those in the G2 group (Figure 6E). Furthermore, the subtypes were significantly related to immune checkpoints, with the G1 group having a greater abundance of immune checkpoints than the G2 group (Figure 6G). The chemotherapy response in the two subtypes was further explored by comparing the IC₅₀s of doxorubicin and cisplatin. A significant decrease in chemotherapy response was observed in G1 compared with G2 based on IC50 of doxorubicin and cisplatin (Figures 6F, H).

FIGURE 6

Figure 6 The subtypes had a correlation with immune and immune response. (A) A comparison of LRG expression in the two subtypes was also carried out. (B) With a better clinical prognosis, the G2 group had lower mRNAsi scores. (C) LRGs were associated with immune infiltration by the EPIC algorithm. (D) TIDE scores of group G1 were significantly higher than those of group G2. (E) The subtypes had a significant association with immune cell. (G) The G1 group having a greater abundance of immune checkpoints than the G2 group. (F, H) The chemotherapy response in the two subtypes was further explored by comparing the IC₅₀s of doxorubicin and cisplatin. Statistical analyses were conducted by t-test, Pearson correlation, log-rank test, and ANOVA. p < 0.05 was considered statistically significant. *p < 0.05; **p < 0.01; ***p < 0.001.

3.6 HPS4 as a potential therapeutic target in LIHC

A univariate Cox regression model was used to develop the nomogram (Figure 7A). On the basis of its c-index of 0.71, it displayed relatively good predictive ability (Figure 7B). According to the calibration plots, the predicted OS and observed OS at 1, 3, and 5 years were well concordant (Figure 7C). As HPS4 and a worse prognosis are strongly associated with LIHC, we further explored whether HPS4 contributes to liver cancer proliferation. SK-hep-1 and Huh7 cells were treated with si-HPS4 to interfere with HPS4 expression (Figure 7D). Moreover, apoptosis could contribute to tumor cell growth; the level of apoptosis in liver cancer cells was also detected using flow cytometry. A significant increase in apoptotic cells was observed when HPS4 was knocked down by si-HPS4 (Figure 7E). The interference of HPS4 with SK-hep-1 and Huh7 cells significantly decreased proliferation by the colony formation assay (Figures 7F). Furthermore, we analyzed the association between HPS4 and a proliferation marker (MKI67). We found that HPS4 was positively related to MKI67 (Figures 7G). Collectively, HPS4 may be a promising therapeutic target and biomarker for LIHC, and regulating its activity may be an effective treatment strategy.

FIGURE 7

Figure 7 Identifying HPS4 as novel therapy target for LIHC. (A) A univariate Cox regression model was used to develop the nomogram. (B) The nomogram displayed relatively good predictive ability. (C) The predicted OS and observed OS at 1, 3, and 5 years. (D) SK-hep-1 and Huh7 cells were treated with si-HPS4 to interfere with HPS4 expression. (E) A significant increase in apoptotic cells was observed when HPS4 was knocked down by si-HPS4. (F) The interference of HPS4 with SK-hep-1 and Huh7 cells significantly decreased proliferation by the colony formation assay. (G) The association between HPS4 and a proliferation marker (MKI67). Statistical analyses included Cox regression, log-rank test, t-test, and Pearson correlation. p < 0.05 was considered statistically significant. ***p < 0.001.

3.7 Analysis of single-cell RNA sequencing data

We further explored the relation between HPS4 and immune cells based on eight single-cell datasets. Our results demonstrated that HPS4 had a significant correlation with CD8 T cell, B cell, and macrophage (Figure 8A). We next analyzed the expression of HPS4 in the hepatocyte population. We found that HPS4 was expressed in cluster 1, cluster 3, cluster 6, cluster 8, cluster 10, cluster 11, and cluster 17. Their corresponding cell populations are activated T cell, sinusoidal endothelial cell, myeloid cell, dendritic cell, smooth muscle cell, hepatocyte_APOA2 high, and epithelial cell_SCGB3A1 high (Figures 8B–D). Similarly, we found that HPS4 was also expressed in hepatocellular carcinoma cells including HEP3B217, HUH6, JHH6, JHH7, LI7, SNU423, and SNU449 (Figures 8E–G). The gene expression level was explored in all cancer single-cell samples. HPS4 was highly expressed in GBM, LUAD, PDAC, and AML relative to other tumors (Figure 8H).

FIGURE 8

Figure 8 The analysis of single-cell RNA sequencing. (A) HPS4 had a significant correlation with CD8 T, B cell, and macrophage based on eight single-cell datasets. (B–D) HPS4 was expressed in cluster 1, cluster 3, cluster 6, cluster 8, cluster 10, cluster 11, and cluster 17. Their corresponding cell populations are activated T cell, sinusoidal endothelial cell, myeloid cell, dendritic cell, smooth muscle cell, hepatocyte_APOA2 high, and epithelial cell_SCGB3A1 high. (E–G) HPS4 was also expressed in hepatocellular carcinoma cells including HEP3B217, HUH6, JHH6, JHH7, LI7, SNU423, and SNU449. (H) The HPS4 expression level was explored in all cancer single-cell samples. Correlation, clustering, and differential expression analyses were performed. p < 0.05 was considered statistically significant.

4 Discussion

Globally, LIHC ranks sixth among the most common cancers. Despite the application of new treatment strategies for LIHC, efficacy remains unsatisfactory. Thus, identification of specific prognostic markers is essential for guiding LIHC therapy and improving OS of patients. It has been reported that lysosome functions are drastically altered during cancer progression, including alterations in lysosome volume, localization, and composition within the cell. A number of cancers show significant overexpression of lysosomal hydrolases, which increases invasion of the tumor (11). It is therefore imperative to understand how LRGs affect LIHC and determine if they serve as prognostic indicators. There has not been any reporting on a prognostic model based on LRGs.

To verify the efficacy of lysosome-related genes on the prognosis of hepatocellular carcinoma, we identified 14 genes out of 61 lysosome-related genes (LRGs) as significant predictors in LIHC. The 14 genes are ANKRD27, AP3M1, BCL10, CD63, CTSC, GLA, HPS1, HPS4, NPC2, PPT1, RAB7A, TPP1, VPS45, and RAMP3. The 14 lysosome-related genes were screened to classify molecular subgroups by NMF consensus clustering. The patients of LIHC from TCGA data were divided into C1 and C2 clusters. We found that cluster 2 had a worse clinical prognosis compared with cluster 1. In addition, we explored that most common CNVs are heterozygous amplifications and deletions of genes based on the LRGs. The impact of tumor CNV levels on immune evasion is particularly intriguing. We also analyzed the differences in immune infiltration between gene-set CNV groups. Further analysis of the pathway prediction and methylation was performed. An OS prognostic signature was constructed using LASSO regression analysis based on the LRGs. As an external verification database, we used the ICGC-LIHC cohorts for training the prognosis model.

Tumor-associated immune responses are regulated and determined by immune cells in the tumor microenvironment (25–27). We also explored the association between LRGs and immune cells. According to some studies, declining CD4+ T cell counts can contribute to liver cancer development (28, 29), whereas our results demonstrated that ANKRD27, AP3M1, BCL10, CTSC, HPS1, HPS4, PPT1, RAB7A, TPP1, and RAMP3 had a negative relation to immune infiltration of CD4 T cells using the EPIC algorithm. We found that subtypes had a significant association with immune cells including B cell, CD4 T cell, CD8 T cell, endothelial cell, macrophage, and NK cell. Several cancer types show a correlation between macrophage density and poor prognosis, suggesting that macrophages play a key role in tumor progression (30–32). Macrophages in the G1 group had lower immune scores than those in the G2 group. Furthermore, the TIDE algorithm showed low efficacy of immune checkpoint blockade therapy (ICB) based on high TIDE scores (21, 33). As a result, the TIDE scores of group G1 were significantly higher than those of group G2. A profound advance in cancer therapy has come from the discovery that immunocheckpoint molecules overexpress in the tumor microenvironment and contribute to antitumor immunity evasion (34–36). Furthermore, the subtypes were significantly related to immune checkpoints, with the G1 group having a greater abundance of immune checkpoints than the G2 group. The resistance to chemotherapy contributes significantly to cancer mortality (37, 38). Thus, the chemotherapy response in the two subtypes was further explored by comparing the IC₅₀s of doxorubicin and cisplatin. A significant decrease in chemotherapy response was observed in G1 compared with G2 based on the IC₅₀ of doxorubicin and cisplatin.

We developed the nomogram by the univariate Cox regression model. Our result demonstrated that HPS4 and a worse prognosis are strongly associated with LIHC. SK-hep-1 and Huh7 cells were treated with si-HPS4 to interfere with HPS4 expression. The interference of HPS4 with SK-hep-1 and Huh7 cells significantly decreased proliferation and colony formation. Moreover, apoptosis could contribute to tumor cell growth; the level of apoptosis in liver cancer cells was also detected using flow cytometry. A significant increase in apoptotic cells was observed when HPS4 was knocked down by si-HPS4. We found that HPS4 was positively related to MKI67 whereas MKI67 is an important marker of proliferation (39). As a tool for studying cell types and states, single-cell RNA sequencing (scRNA-seq) is becoming increasingly important. We further explored the relation between HPS4 and immune cells based on eight single-cell datasets. Our results demonstrated that HPS4 had a significant correlation with CD8 T, B cell, and macrophage. We found that the genes were mainly expressed in activated T cell, sinusoidal endothelial cell, myeloid cell, dendritic cell, smooth muscle cell, hepatocyte_APOA2 high, and epithelial cell_SCGB3A1 high. Collectively, it turns out that HPS4 may be a promising therapeutic target and biomarker for LIHC, and regulating its activity may be an effective treatment strategy.

The study has some limitations. There is still a lack of understanding of how HPS4 impacts proliferation or apoptosis of liver cancer cells. This is our next endeavor to further explore how HPS4 regulates LIHC in vivo and in vitro.

5 Conclusion

A prognostic model based on LRGs was developed in this study to predict the prognosis of patients with LIHC. HPS4, which could affect proliferation and apoptosis of liver cancer cells, may be a promising therapeutic target and biomarker for LIHC. HPS4 may be a promising therapeutic target and biomarker for LIHC, and regulating its activity may be an effective treatment strategy.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

K-JH designed this work. K-JH performed bioinformatic analyses. K-JH wrote the manuscript. K-JH and ZN revised the manuscript and supervised the whole experiment. All authors have read and approved the final submitted manuscript.

Funding

This project was funded by Quzhou Municipal Science and Technology Bureau (2023K117) and Quzhou Medical and Health personnel scientific research fund project (KYQD2022-18).

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2023.1221498/full#supplementary-material

Supplementary Figure 1 | LRGs were associated with drug sensitivity based on the CTRP (a) and GDSC resource (b). Statistical analysis was performed using Pearson correlation. P < 0.05 was considered statistically significant. *P < 0.05; **P < 0.01; *** P < 0.001.

Supplementary Figure 2 | The association between subtypes and ferroptosis and m6A. Enrichment analysis and correlation tests were conducted. P < 0.05 was considered statistically significant. *P < 0.05; **P < 0.01; *** P < 0.001.

Supplementary Figure 3 | The differential genes between cluster G1 and G2; (a) The volcano plot; (b) heat map of differential genes; (c) Pathway enrichment analysis. P < 0.05 was considered statistically significant. *P < 0.05; **P < 0.01; *** P < 0.001.

Supplementary Figure 4 | Correlation between high and low risk groups and immune cells. Correlation analysis was conducted using Pearson correlation. P < 0.05 was considered statistically significant. *P < 0.05; **P < 0.01; *** P < 0.001.

References

1. Akinyemiju T, Abera S, Ahmed M, Alam N, Alemayohu MA, Allen C, et al. The burden of primary liver cancer and underlying etiologies from 1990 to 2015 at the global, regional, and national level: results from the global burden of disease study 2015. JAMA Oncol (2017) 3(12):1683–91.

PubMed Abstract | Google Scholar

2. Valery PC, Laversanne M, Clark PJ, Petrick JL, McGlynn KA, Bray F, et al. Projections of primary liver cancer to 2030 in 30 countries worldwide. Hepatology (2018) 67(2):600–11. doi: 10.1002/hep.29498

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Tseng HC, Xiong W, Badeti S, Yang Y, Ma M, Liu T, et al. Efficacy of anti-CD147 chimeric antigen receptors targeting hepatocellular carcinoma. Nat Commun (2020) 11(1):4810. doi: 10.1038/s41467-020-18444-2

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Marquardt JU, Andersen JB, Thorgeirsson SS. Functional and genetic deconstruction of the cellular origin in liver cancer. Nat Rev Cancer (2015) 15(11):653–67. doi: 10.1038/nrc4017

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Yan QH, Xu DG, Shen YF, Yuan DL, Bao JH, Li HB, et al. Observation of the effect of targeted therapy of 64-slice spiral CT combined with cryoablation for liver cancer. World J Gastroenterol (2017) 23(22):4080–9. doi: 10.3748/wjg.v23.i22.4080

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Tao J, Zhang R, Singh S, Poddar M, Xu E, Oertel M, et al. Targeting β-catenin in hepatocellular cancers induced by coexpression of mutant β-catenin and K-Ras in mice. Hepatology (2017) 65(5):1581–99. doi: 10.1002/hep.28975

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wischhusen JC, Chowdhury SM, Lee T, Wang H, Bachawal S, Devulapally R, et al. Ultrasound-mediated delivery of miRNA-122 and anti-miRNA-21 therapeutically immunomodulates murine hepatocellular carcinoma in vivo. J Control Release (2020) 321:272–84. doi: 10.1016/j.jconrel.2020.01.051

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Wu R, Murali R, Kabe Y, French SW, Chiang YM, Liu S, et al. Baicalein targets GTPase-mediated autophagy to eliminate liver tumor-initiating stem cell-like cells resistant to mTORC1 inhibition. Hepatology (2018) 68(5):1726–40. doi: 10.1002/hep.30071

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Barraud L, Merle P, Soma E, Lefrançois L, Guerret S, Chevallier M, et al. Increase of doxorubicin sensitivity by doxorubicin-loading into nanoparticles for hepatocellular carcinoma cells in vitro and in vivo. J Hepatol (2005) 42(5):736–43. doi: 10.1016/j.jhep.2004.12.035

PubMed Abstract | CrossRef Full Text | Google Scholar

10. OuYang HY, Xu J, Luo J, Zou RH, Chen K, Le Y, et al. MEP1A contributes to tumor progression and predicts poor clinical outcome in human hepatocellular carcinoma. Hepatology (2016) 63(4):1227–39. doi: 10.1002/hep.28397

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kundu ST, Grzeskowiak CL, Fradette JJ, Gibson LA, Rodriguez LB, Creighton CJ, et al. TMEM106B drives lung cancer metastasis by inducing TFEB-dependent lysosome synthesis and secretion of cathepsins. Nat Commun (2018) 9(1):2731. doi: 10.1038/s41467-018-05013-x

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hui Y, Yi X, Wibowo D, Yang G, Middelberg APJ, Gao H, et al. Nanoparticle elasticity regulates phagocytosis and cancer cell uptake. Sci Adv (2020) 6(16):eaaz4316. doi: 10.1126/sciadv.aaz4316

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zhang X, Anthony B, Chai Z, Lee Dobbins A, Sutton RB, Li C, et al. Membrane fusion FerA domains enhance adeno-associated virus vector transduction. Biomaterials (2020) 241:119906. doi: 10.1016/j.biomaterials.2020.119906

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Bonam SR, Wang F, Muller S. Lysosomes as a therapeutic target. Nat Rev Drug Discovery (2019) 18(12):923–48. doi: 10.1038/s41573-019-0036-1

CrossRef Full Text | Google Scholar

15. Lee CH, Song DK, Park CB, Choi J, Kang GM, Shin SH, et al. Primary cilia mediate early life programming of adiposity through lysosomal regulation in the developing mouse hypothalamus. Nat Commun (2020) 11(1):5772. doi: 10.1038/s41467-020-19638-4

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Berg RD, Levitte S, O'Sullivan MP, O'Leary SM, Cambier CJ, Cameron J, et al. Lysosomal disorders drive susceptibility to tuberculosis by compromising macrophage migration. Cell (2016) 165(1):139–52. doi: 10.1016/j.cell.2016.02.034

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhang L, Gao Y, Zhang X, Guo M, Yang J, Cui H, et al. TSTA3 facilitates esophageal squamous cell carcinoma progression through regulating fucosylation of LAMP2 and ERBB2. Theranostics (2020) 10(24):11339–58. doi: 10.7150/thno.48225

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Maehr R, Mintern JD, Herman AE, Lennon-Duménil AM, Mathis D, Benoist C, et al. Cathepsin L is essential for onset of autoimmune diabetes in NOD mice. J Clin Invest (2005) 115(10):2934–43. doi: 10.1172/JCI25485

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Khalkhali-Ellis Z, Hendrix MJ. Elucidating the function of secreted maspin: inhibiting cathepsin D-mediated matrix degradation. Cancer Res (2007) 67(8):3535–9. doi: 10.1158/0008-5472.CAN-06-4767

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Runkle KB, Meyerkord CL, Desai NV, Takahashi Y, Wang HG, et al. Bif-1 suppresses breast cancer cell migration by promoting EGFR endocytic degradation. Cancer Biol Ther (2012) 13(10):956–66. doi: 10.4161/cbt.20951

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res (2021) 49(D1):D1420–d1430. doi: 10.1093/nar/gkaa1020

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Han X, Zhou Z, Fei L, Sun H, Wang R, Chen Y, et al. Construction of a human cell landscape at single-cell level. Nature (2020) 581(7808):303–9. doi: 10.1038/s41586-020-2157-4

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Database resources of the national genomics data center, China national center for bioinformation in 2022. Nucleic Acids Res (2022) 50(D1):D27–d38.

PubMed Abstract | Google Scholar

25. Nagarsheth N, Wicha MS, Zou W. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol (2017) 17(9):559–72. doi: 10.1038/nri.2017.49

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Riscal R, Skuli N, Simon MC. Even cancer cells watch their cholesterol! Mol Cell (2019) 76(2):220–31. doi: 10.1016/j.molcel.2019.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Chen X, Litzenburger UM, Wei Y, Schep AN, LaGory EL, Choudhry H, et al. Joint single-cell DNA accessibility and protein epitope profiling reveals environmental regulation of epigenomic heterogeneity. Nat Commun (2018) 9(1):4590. doi: 10.1038/s41467-018-07115-y

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Dal Maso L, Polesel J, Serraino D, Lise M, Piselli P, Falcini F, et al. Pattern of cancer risk in persons with AIDS in Italy in the HAART era. Br J Cancer (2009) 100(5):840–7. doi: 10.1038/sj.bjc.6604923

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Weng JJ, Wei JZ, Li M, Lu JL, Qin YD, Jiang H, et al. Effects of hepatitis B virus infection and antiviral therapy on the clinical prognosis of nasopharyngeal carcinoma. Cancer Med (2020) 9(2):541–51. doi: 10.1002/cam4.2715

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Terashima Y, Toda E, Itakura M, Otsuji M, Yoshinaga S, Okumura K, et al. Targeting FROUNT with disulfiram suppresses macrophage accumulation and its tumor-promoting properties. Nat Commun (2020) 11(1):609. doi: 10.1038/s41467-020-14338-5

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Guerin MV, Regnier F, Feuillet V, Vimeux L, Weiss JM, Bismuth G, et al. TGFβ blocks IFNα/β release and tumor rejection in spontaneous mammary tumors. Nat Commun (2019) 10(1):4131.

PubMed Abstract | Google Scholar

32. LaGory EL, Giaccia AJ. The ever-expanding role of HIF in tumour and stromal biology. Nat Cell Biol (2016) 18(4):356–65. doi: 10.1038/ncb3330

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Kim MH, Kim JH, Lee JM, Choi JW, Jung D, Cho H, et al. Molecular subtypes of oropharyngeal cancer show distinct immune microenvironment related with immune checkpoint blockade response. Br J Cancer (2020) 122(11):1649–60. doi: 10.1038/s41416-020-0796-8

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Altan M, Pelekanou V, Schalper KA, Toki M, Gaule P, et al. B7-H3 expression in NSCLC and its association with B7-H4, PD-L1 and tumor-infiltrating lymphocytes. Clin Cancer Res (2017) 23(17):5202–9. doi: 10.1158/1078-0432.CCR-16-3107

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Khasraw M, Reardon DA, Weller M, Sampson JH. PD-1 Inhibitors: Do they have a Future in the Treatment of Glioblastoma? Clin Cancer Res (2020) 26(20):5287–96. doi: 10.1158/1078-0432.CCR-20-1135

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Barroso-Sousa R, Barry WT, Garrido-Castro AC, Hodi FS, Min L, Krop IE, et al. Incidence of endocrine dysfunction following the use of different immune checkpoint inhibitor regimens: A systematic review and meta-analysis. JAMA Oncol (2018) 4(2):173–82. doi: 10.1001/jamaoncol.2017.3064

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bester AC, Lee JD, Chavez A, Lee YR, Nachmani D, Vora S, et al. An integrated genome-wide CRISPRa approach to functionalize lncRNAs in drug resistance. Cell (2018) 173(3):649–664.e20. doi: 10.1016/j.cell.2018.03.052

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ohkuma R, Yada E, Ishikawa S, Komura D, Kubota Y, Hamada K, et al. High levels of human epididymis protein 4 mRNA and protein expression are associated with chemoresistance and a poor prognosis in pancreatic cancer. Int J Oncol (2021) 58(1):57–69.

PubMed Abstract | Google Scholar

39. Fougner C, Bergholtz H, Norum JH, Sørlie T. Re-definition of claudin-low as a breast cancer phenotype. Nat Commun (2020) 11(1):1787. doi: 10.1038/s41467-020-15574-5

PubMed Abstract | CrossRef Full Text | Google Scholar