Predicting the prognosis of hepatocellular carcinoma based on genes related to polyamine metabolism

Data collection and acquisition

In this study, the bulk dataset (training set) containing the RNA-seq data of The Cancer Genome Atlas (TCGA)-HCC was downloaded through TCGA Genomic Data Commons (GDC) Application Programming Interface (API) (https://portal.gdc.cancer.gov/). Samples without clinical follow-up data or clearly defined statuses were eliminated. The Ensembl IDs were transformed to gene symbols, while the average expression of genes with multiple gene symbols was calculated. After screening, a total of 370 primary tumor samples and 50 adjacent non-tumor samples were successfully acquired. The International Cancer Genome Consortium Liver Cancer Japan Project (ICGC-LIRI-JP) dataset, a bulk dataset obtained from A Database of Hepatocellular Carcinoma Expression Atlas (HCCDB, http://lifeome.net/database/hccdb/), included 212 liver cancer samples and used as a validation set.

GSE166635 containing the scRNA-seq data of two HCC tumor samples (sample numbers: GSM5076749 and GSM5076750) was downloaded from the Gene Expression Omnibus (GEO) database. A total of 59 PMRGs were collected from the “REACTOME_METABOLISM_OF_POLYAMINES” pathway in the MSigDB database (http://www.gsea-msigdb.org/).

The scRNA-seq data preprocessing

The scRNA-seq dataset GSE166635 was preprocessed under rigorous screening criteria. Specifically, each gene should be expressed in three cells at least, and each individual cell expressed a minimum of 200 genes. Cells were retained if they had 200 to 8,000 genes (nFeature_RNA), <100,000 transcripts (nCount_RNA), and <15% mitochondrial gene expression (percent.mito). Then, the NormalizeData function (Stuart et al., 2019) was used to perform log-transformation on the data, and the FindVariableFeatures function (Stuart et al., 2019) was applied to accurately filter highly variable genes. The expression values of all genes were normalized by the ScaleData function (Stuart et al., 2019). Utilizing the RunPCA function, principal component analysis (PCA) was performed (Stuart et al., 2019). Batch effects among various samples were removed applying the harmony package (Korsunsky et al., 2019), selecting the top 20 PCs for dimensionality reduction by UMAP. Then the cells were clustered into subsets by the FindNeighbors and FindClusters (Stuart et al., 2019) and finally annotated to specific cell types utilizing the marker genes provided by the CellMarker2.0 database (Hu et al., 2022).

Single-sample GSEA

The enrichment scores of PMRGs in the TCGA-HCC cohort were computed by single-sample GSEA (ssGSEA) (Barbie et al., 2009). Specifically, the relative enrichment of this gene set in each sample was estimated by comparing the gene expression data of each sample in the TCGA-HCC cohort to a specific set of PMRGs. The enrichment score calculated by ssGSEA indicated the extent to which a particular gene set within a sample was either synchronously upregulated or downregulated.

Validation of molecular subtypes classified using the PMRGs

Univariate Cox regression analysis was utilized to elucidate the relationship between PMRGs and the prognosis of HCC. By setting the criterion of p < 0.05, PMRGs with significant prognostic significance were screened to classify molecular subtypes of HCC. Subsequently, the ConsensusClusterPlus method (Wilkerson & Hayes, 2010) was employed to construct a consensus matrix for accurate molecular subtyping. The clustering process utilized the “km” algorithm and 500 bootstraps, using “1-Spearman correlation” as the distance metric with 80% of the training set patients in each bootstrap. The cluster number was between 2 and 10 for analysis. Prognostic differences among various molecular subtypes were compared by the Kaplan–Meier (KM) curves, and differences in clinicopathological features (age (>60/≤60), stage (I–IV), gender (male/female), T stage (T1–T4), and grade (G1–G4)) of different molecular subtypes in the TCGA-HCC dataset were also explored.

Discovery of DEGs by enrichment analyses

The DEGs between the identified subtypes were selected by the limma package (Ritchie et al., 2015; Wang et al., 2024) under the criteria of FDR < 0.05 and |log2FC|>log2(2). Subsequently, using the R package clusterProfiler (Yu et al., 2012), KEGG and GO (in biological process, GO-BP) functional enrichment analyses were performed on these DEGs (Wang et al., 2025).

Biomarkers screening through machine learning

Prognostically significant genes (p < 0.05) for HCC were selected from the DEGs by univariate Cox regression analysis and then subjected to Lasso regression analysis using the glmnet R package (Friedman et al., 2021). Based on the results of Lasso analysis, stepwise multiple regression analysis was further performed. Akaike information criterion (AIC) was utilized to optimize the model during the process. Specifically, the stepAIC method from the MASS package (Ripley et al., 2013) was adopted, iteratively removing one variable at a time beginning with the most complex model while monitoring the changes in AIC value.

Establishment of the risk model

The formula RiskScore = Σβi × Expi was used to compute RiskScore for each patient, with Expi and β presenting the gene expression and the Cox regression coefficient of the gene. Subsequently, the RiskScore was normalized by z-score, and high-risk and low-risk groups were divided by the threshold value of “0”. The KM survival curve was used to conduct prognostic analysis, and significant differences were analyzed by the log-rank test. Furthermore, the timeROC R package (Blanche, Dartigues & Jacqmin-Gadda, 2013) was employed as a prognostic classification tool to conduct a receiver operating characteristic (ROC) analysis for the RiskScore.

Construction of a nomogram

Univariate and multivariate Cox regression analyses were carried out in the TCGA-HCC dataset to explore whether the RiskScore was independent of Stage, T.stage, N.stage, M.stage, Grade, Risktype, gender. Additionally, survival and risk evaluation for patients in this dataset was further improved by developing a nomogram.

Relation between the RiskScore and immune microenvironment

Tumor Immune Estimation Resource (TIMER) (Li et al., 2020), MCP-counter (Becht et al., 2016), and Cell Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) (Chen et al., 2018) were employed to examine the status of various immune cell types, so as to clarify the inherent connection between the RiskScore and immune microenvironment of different patients.

Correlation analysis between drug sensitivity and the RiskScore

The IC₅₀ of drugs for the TCGA-HCC dataset samples was predicted by the oncoPredict R package (Maeser, Gruener & Huang, 2021). Next, the relationship between the RiskScore and drug sensitivity was examined applying by Spearman’s correlation analysis, with a p-value < 0.05 and |cor|>0.4 denoting a significant relationship.

Cell culture

Human liver immortalized cells THLE-2 (C5664), human HCC cells HuH-7 (C5176) were obtained from BDBio (https://www.biocode.cn/) (Hangzhou, China). All the cells used were identified by short tandem repeat (STR) to confirm whether they were contamination-free. Cells were cultured using DMEM (BDBio, Hangzhou, China) with 10% fetal bovine serum (FBS) (C0226; Beyotime, Shanghai, China) and 1% penicillin-streptomycin (P/S) (15140148; Thermo Fisher Scientific, Waltham, MA, USA) and incubated at 37 °C in the concentration of 5% CO₂ and saturated humidity.

RNA extraction and qRT-PCR experiment

Following the instructions, total RNA was separated from THLE-2 and HuH-7 cell samples using Total RNA Extraction Kit (R1200; Solarbio, Beijing, China). RNA quality was assessed by measuring the OD260/OD280 ratio with a NanoDrop microspectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Next, the cDNA template was obtained using BeyoRT™ II M-MLV reverse transcriptase (D7160; Beyotime, Shanghai, China). Changes in the expressions of the biomarkers obtained from computational analysis of THLE-2 and HuH-7 cells were analyzed employing BeyoFast™ SYBR Green qPCR Mix (D7265; Beyotime, Shanghai, China) and qRT-PCR technology on an ABI7300 Realtime PCR System (Thermo Fisher Scientific, Waltham, MA, USA). Gene expression was quantified by 2^−ΔΔCT method and normalized to that of GAPDH. The primers were designed by Origene (Rockville, MD, USA). See Table S1 for the relevant sequences and primers.

Cell transfection

HuH-7 cells were transfected with siRNA by the use of Lipofectamine 2000 Transfection Kit (11668027; Invitrogen, Carlsbad, CA, USA). Following the instructions, Lipofectamine 2000 and siRNA were diluted to 100 nM and 3 μL/well, respectively, to ensure an optimal transfection. Next, the cells were incubated at 37 °C for 48 hours (h). To prevent off-target effects, two siRNA targeting sequences and a negative control (designed and purchased from Merck KGaA, Darmstadt, Germany) with the optimal transfection efficiencies were used for subsequent experiment, which were labeled as si-G6PD#1 (target sequence: 5′-GCCTTCCATCAGTCGGATACA-3′), si-G6PD#2 (target sequence: 5′-GCTGACATCCGCAAACAGAGT-3′), and si-NC.

Cell viability assay

Cell proliferation of HuH-7 cells transfected with si-NC and si-G6PD was detected using CCK-8 assay kit (C0037; Beyotime, Shanghai, China). The cells were inoculated at a density of 5 × 10³ cells/well and 100 µL of medium was added to each well to ensure cell growth in a suitable environment. During cell culture, CCK-8 solution (10 µL) was supplemented to each 96-well plate at different time points (24 h, 48 h, 72 h) and incubated for 30 minutes (min) to react with the cells. Finally, the absorbance was measured at 450 nm.

Cell migration assay

HuH-7 cells were evenly seeded into 6-well plates, with 2 mL of complete medium in each well. The cells were incubated at 37 °C with 5% CO₂ until a confluent monolayer was formed. When the cells reached approximately 90% confluence, a straight scratch was gently created across the cell monolayer using a sterilized 200 µL pipette tip. The existing medium was then removed, and the wells were washed twice with PBS (C0221A; Beyotime, Shanghai, China) to eliminate cell debris. Then fresh serum-free medium was added, and the plates were returned to the incubator under the original conditions to continue the incubation. The plates were removed at 0 and 48 h. Photographs of the scratched area were taken with a microscope (Olympus Corporation, Tokyo, Japan), and the time and abnormality were recorded. Finally, changes in scratch width were measured by ImageJ software and the scratch closure (%) was calculated at each time point by the formula: initial scratch width-current scratch width/initial scratch width × 100%.

Cell invasion capacity

Cell suspension was seeded into the upper Transwell chamber (8 μm, Corning, Corning, NY, USA) at a density of 3 × 10⁴ cells/chamber, while complete medium was added to the lower chamber as a chemoattractant. After 48-h incubation at 37 °C with 5% CO₂, non-invaded cells were removed. The membrane was then washed with PBS (C0221A; Beyotime, Shanghai, China), immobilized in 4% paraformaldehyde (P0099; Beyotime, Shanghai, China) for 15 min, and subsequently colored by crystal violet solution (C0121; Beyotime, Shanghai, China) for 10 min. After rinsing in distilled water, the membrane was air-dried and the invaded cells were quantified from at least five random fields using an inverted microscope (Olympus Corporation, Tokyo, Japan).

Statistical analysis

R software (version 3.6.0) and GraphPad Prism (version 8.0.2) were used for statistical analyses. Student’s t-test or Wilcoxon rank-sum test was used to compare differences in continuous variables between two groups. For datasets involving three or more variables, ANOVA was employed to assess the overall differences among the groups. Subsequently, pairwise comparisons were performed using Sidak’s multiple comparison test. Statistically significant was defined p < 0.05, ns meant not significant.

Molecular subtyping using prognostically significant PMRGs

In the TCGA-HCC cohort, ssGSEA calculation showed that the enrichment scores of PMRGs were higher in the cancer tissue scores than in normal tissues (Fig. 1A, p < 0.05). A total of 42 PMRGs (p < 0.05) were found to be closely associated with the prognostic outcomes of HCC, and most of them were risk genes. Next, based on the expression profiles of these 42 PMRGs, the consensus clustering analysis revealed a comparatively stable clustering result at k = 2, which classified two molecular subtypes (C1 and C2) (Fig. 1B). Notably, the prognosis of C2 was more unfavorable than C1 (Fig. 1C, p < 0.05). Meanwhile, the distribution of various clinical features across different molecular subtypes in the TCGA-HCC dataset was compared. It was found that subtype C1 had a relatively more patient at higher clinical grades (Fig. 1D).

Functional analysis between the two molecular subtypes

We screened 100 DEGs between the subtypes C1 and C2 (Fig. 2A) and analyzed the distribution of these DEGs between the two subtypes (Fig. 2B). KEGG analysis demonstrated that the DEGs were mainly enriched in pathways, including primary bile acid (BA) biosynthesis, PPAR signaling pathway, bile secretion, IL-17 signaling pathway, phenylalanine metabolism (Fig. 2C). GO-BP enrichment analysis showed that these genes were mainly enriched in pathways, including steroid catabolic process, steroid metabolic process, steroid biosynthetic process, cellular hormone metabolic process, BA biosynthetic process (Fig. 2D).

Development of a clinical prognostic model and verification

Univariate COX regression analysis identified 70 genes that may markedly influence HCC prognosis (p < 0.05). To optimize the model, Lasso regression demonstrated that variable coefficients gradually approached zero with increasing λ values (Fig. 3A). The optimal performance of the model was achieved at lambda = 0.0674 with 10-fold cross-validation, identifying six key genes (Fig. 3B). Further stepwise multiple regression analysis reduced the six genes to four genes (Fig. 3C), namely, G6PD, S100A9, AKR1B15, and ADH4. A risk model was established as: RiskScore = 0.188*G6PD + 0.093*AKR1B15 + 0.079*S100A9 − 0.059*ADH4.

Subsequently, each sample in the training set was assigned with a RiskScore, normalized by z-score and classified into high- and low-risk groups, which contained 157 and 213 samples, respectively (Fig. 3D). The ROC curve showed a high AUC value of the RiskScore, indicating its robust predictive performance (Fig. 3E). Next, comparison on the overall survival (OS) rates of the patients with different RiskScores demonstrated a negative correlation between a worse survival and a high RiskScore (Fig. 3F, p < 0.05). In addition, verification in the validation set cohort showed similar results to those obtained in the training set (Figs. 3G, 3H, and 3I).

Results of in vitro experiments

Cellular validation assays demonstrated that the expressions of G6PD, AKR1B15 and S100A9 were all notably elevated in HuH-7 cells (p < 0.05, Fig. 4A), while the level of ADH4 was not significantly different in HuH-7 cells and THLE-2 cells (p > 0.05, Fig. 4A). Among the four biomarkers, G6PD showed the most significant differential expression, and the transfection efficiency of si-G6PD#2 was higher than si-G6PD#1. Thus, si-G6PD#2 was selected for subsequent experiments (p < 0.05, Fig. 4B). CCK8 experiments revealed that the cell viability was notably higher in the si-NC group than in the si-G6PD#2 group at 24 h, 48 h and 72 h (p < 0.05, Fig. 4C). Wound healing and Transwell assays also indicated that the cell migratory and invasive abilities of the si-NC group were higher than those of the si-G6PD#2 group (p < 0.05, Figs. 4D and 4E). These findings suggested that the silencing of G6PD gene expression suppressed the viability, migratory and invasive abilities of HCC cells.

Development of a nomogram

The RiskScore was verified as the most important prognosis factor (p < 0.05) in the TCGA-HCC dataset by univariate and multivariate Cox regression analysis (Figs. 5A and 5B). After developing a nomogram, the RiskScore also manifested the strongest influence on the OS prediction in HCC (Fig. 5C). Further, the calibration curve calibration curves were plotted to test the prediction accuracy of the nomogram, which demonstrated close alignment between the estimated and actual outcomes at 1-, 3-, and 5-year time points (Fig. 5D). This indicated an accurate prediction ability of the nomogram. In addition, the model reliability was assessed by decision curve analysis (DCA). We found that the nomogram also showed the strongest capability in estimating the OS in HCC than several other clinicopathological features (Fig. 5E).

Analysis of drug sensitivity and immune microenvironment

The association between immune microenvironment and the RiskScore was comprehensively examined. The results of TIMER indicated that the infiltration of neutrophil cells, dendritic cells (DCs), T cell CD8, macrophages, B cells, T cells CD4 in the high-risk group was significantly higher than those in the control group (p < 0.05, Fig. 6A). The MCP-counter further indicated markedly higher scores of most immune cells such as monocytic lineage, T cells, myeloid DCs, CD8 T cells, fibroblasts, B lineages, cytotoxic lymphocytes in the high-risk group than those in the low-risk group (p < 0.05, Fig. 6B). CIBERSORT analysis showed that monocytes, B cells naïve, macrophages M1, T cells CD4 memory resting, mast cells resting, T cells gamma delta, etc. had markedly higher scores of in the low-risk group, whereas the high-risk group had markedly higher scores of T cell regulatory Tregs, macrophage M0, T cell follicular helper, DC resting, etc. (p < 0.05, Fig. 6C).

Subsequently, the IC₅₀ values of each drug in the samples of the TCGA-HCC dataset were calculated to identify drugs with significant differences. A significant correlation was found between 13 drugs and the RiskScore. In particular, SB505124_1194 and Doramapimod_1042 were positively correlated with the RiskScore, while the other 11 drugs such as lapatinib_1558 were all negatively correlated with the RiskScore (Fig. 6D).

Singl-cell atlas

To elucidate the expression distribution of the key genes in different HCC cell types, the scRNA-seq data from two HCC samples (GSM5076749 and GSM5076750) in GSE166635 were analyzed (Figure S1). After cell filtering, normalization, dimensionality reduction, and clustering, a total of 18,369 cells were retained and we clustered 10 major cell subsets. According to the expressions of marker genes in the clusters, eight cell types were annotated, namely T cells, macrophages, cancer-associated fibroblasts (CAFs), B cells, plasma cells, thymic epithelial cells (TECs), mast cells, and HPCs (Figs. 7A and 7B). These 8 cell types were confirmed by identifying their specific marker genes (Fig. 7C). The proportions of cells in the two HCC samples (GSM5076749 and GSM5076750) were different, but both the two were dominated by T cells and macrophages (Fig. 7D). In addition, by analyzing the expressions of the four genes in different cell gene clusters, it was found that both S100A9 and G6PD genes were high-expressed in HPCs and macrophages (Fig. 7E).