Construction of a prognostic model based on protein post-translational modification genes for prediction of immune characteristics and therapeutic response in hepatocellular carcinoma.

Introduction
Hepatocellular carcinoma (HCC) is a leading malignancy of the digestive system, accounting for a significant portion of global cancer mortality [1]. The primary etiologies of HCC include chronic viral hepatitis, cirrhosis, and non-alcoholic steatohepatitis [2]. Due to the asymptomatic nature of early-stage HCC, most patients are diagnosed at an advanced stage, where metastasis has already occurred, leading to a poor prognosis. Current therapeutic strategies, such as hepatic resection, liver transplantation, and ablation, have limited efficacy in advanced cases and are associated with high recurrence rates [3, 4]. Therefore, identifying novel biomarkers that can accurately predict patient prognosis and response to treatment is crucial [5–8].
Immune infiltration plays a key role in cancer progression and patient outcomes. Immune cells, such as T cells, macrophages, and natural killer cells, mediate an anti-tumor response, but tumor cells develop mechanisms to evade immune surveillance [9]. Immunotherapy, particularly immune checkpoint inhibitors (e.g., PD-L1 and CTLA-4 inhibitors), has emerged as a promising therapeutic option. However, due to tumor resistance, not all patients respond effectively to these treatments [10]. The inhibitors of PD-L1 and CTLA4 have been widely employed in immune therapy for cancer patients, offering new hope for cancer treatment. However, due to the development of tumor resistance, immune therapy is not universally effective for all patients [4], therefore, the prediction of the response to immunotherapy will be beneficial to the clinical treatment of patients.
Protein post-translational modifications (PTMs) significantly influence protein function and are closely associated with tumorigenesis and immune responses. PTMs such as phosphorylation, acetylation, and glycosylation have been shown to play critical roles in the development of HCC [11]. Recent research indicates a close association between PTMs and the onset and progression of cancer. For example, PGK1 promotes HCC proliferation by phosphorylating PRAS40 in order to inhibit autophagy-mediated cell death [12]. Melatonin inhibits the activity of MMP9 by histone acetylation, thereby suppressing TPA-induced oral cancer migration [13], and GALNT1 promotes the progression of gastric cancer by enhancing CD44 glycosylation to activate the Wnt/β-catenin signaling pathway [14]. In recent years, a growing body of evidence has highlighted the central role of PTMs in reshaping the tumor immune microenvironment, positioning them as key mechanisms for overcoming therapeutic bottlenecks. For example, aerobic glycolysis enhances the expression of PD-L1 through HK2-mediated IκBα phosphorylation, thereby promoting tumor immune evasion [15]. The E3 ubiquitin ligase UBR5 induces immune-suppressing macrophages and thus promoting the growth and metastasis of ovarian cancer [16]. Down-regulation of B3GNT3 inhibits PD-L1 glycosylation and reduces PD-1-PD-L1 interactions, which in turn enhances the killing effect of cytotoxic T cells on triple-negative breast cancer [17]. Furthermore, FOXP3 acetylation modulates immune suppression by regulating Treg cells [18, 19]. USP22 enhances PD-L1 expression through ubiquitination, thereby suppressing the immune microenvironment by depleting infiltrating CD8 + T cells [19, 20]. Thus, PTMs play a crucial role in regulating immune checkpoints like PD-L1. Given that current immunotherapies primarily target these checkpoints, investigating PTMs offers new insights into treatment sensitivity and resistance. Although the importance of PTMs is recognized, there is still a gap in systematically constructing a prognostic model based on PTM gene profiles and comprehensively evaluating its association with immunotherapy response—thereby underscoring the necessity of this study.
This study focuses on seven representative post-translational modifications (PTMs): phosphorylation, acetylation, methylation, SUMOylation, neddylation, ubiquitination, and glycosylation, as their roles in HCC pathogenesis and immune regulation have been well established. Phosphorylation and acetylation are extensively studied modification types in HCC, influencing key signaling pathways such as PI3K/AKT and STAT3, and regulating transcriptional activity through histone modifications [21, 22]. Ubiquitination and SUMOylation are crucial for protein degradation and stability, modulating oncogenes and tumor suppressors (e.g., p62, LKB1) [23, 24]. Glycosylation alters the function of cell surface proteins, including immune checkpoint molecules like PD-L1, thereby influencing immune evasion [25]. Methylation regulates gene expression by controlling methylation and demethylation, subsequently affecting HCC progression. Neddylation is abnormally activated in HCC, altering substrate properties to activate signaling pathways that influence cell proliferation, cycle progression, apoptosis, and other processes, thereby inducing liver disease progression to HCC [26]. Collectively, these seven protein modification types constitute a broad spectrum of regulatory mechanisms highly relevant to HCC biology and therapeutics, making them the foundation for constructing a comprehensive prognostic model. The aim of this study was to construct a prognostic model for HCC based on PTM genes to predict patient prognosis and response to immunotherapy. By collecting 4210 PTM-related genes from public databases and utilizing Weighted Gene Co-expression Network Analysis (WGCNA), key module genes were identified. Using 101 machine learning algorithm combinations, the optimal prognostic model was developed and validated with external datasets. Additionally, the study explored differences in biological functions, immune status, mutation burden, and sensitivity to immunotherapy and chemotherapy between high-risk and low-risk patient groups.

Materials and methods

Data sources and processing
The flowchart for this study is shown in Fig. 1.Transcriptomic and clinical data from the TCGA-LIHC cohort, which includes 50 normal samples and 374HCC samples, were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Additionally, transcriptomic and clinical data for 243 HCC samples from the LIRI-JP cohort were downloaded from the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/). The GSE116174 microarray dataset, consisting of transcriptomic and clinical data for 64 HCC samples, was retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Single-cell transcriptomic data for HCC were sourced from the China National GeneBank (CNSA: CNP0000650), comprising 12 primary HCC samples and 6 recurrent HCC samples. PTM genes were obtained from the MsigDB database (http://software.broadinstitute.org/gsea/index.jsp).

WGCNA analysis
WGCNA was employed to identify key genes associated with HCC. The optimal soft threshold was calculated using the “PickSoftThreshold” function. A co-expression network was constructed using a one-step approach, with the minimum module size set to 50 genes. Pearson correlation was used to assess the relationship between each module and the tumor group, with the module showing the highest correlation designated as the key module.

Differential expression analysis and prognosis analysis of PTM
Differential expression analysis was conducted on the PTM genes within the key module, with differentially expressed genes (DEGs) identified using the criteria |log2FC| > 1.5 and FDR < 0.05. Univariate Cox regression analysis was subsequently applied to the DEGs to select prognostic genes, using a significance threshold of p < 0.001.

Cluster analysis of prognostic PTM genes
Clustering of HCC samples based on the expression of prognostic PTM genes was performed using the “ConsensusClusterPlus” package in R software. Heatmaps of prognostic gene expression across clusters were generated using the “pheatmap” package. Kaplan-Meier survival analysis was conducted to compare survival outcomes between clusters, while the “GSVA” package was used to explore functional differences between clusters.

Construction of prognostic model
We employed ten machine learning algorithms, including least absolute shrinkage and selection operator (LASSO), random forest (RSF), gradient boosting machine (GBM), Cox partial least squares regression (plsRcox), ridge regression, supervised principal component (SuperPC), CoxBoost, survival support vector machine (Survival-SVM), stepwise Cox regression, and elastic net (Enet). These ten machine learning algorithms can be combined into 101 algorithm combinations through either standalone application or pairwise combinations. Among these algorithm combinations, the first algorithm is used to screen for feature genes, while the second algorithm utilizes the screened feature genes to construct a prognostic model. The optimal model combination for constructing the prognostic model was selected based on the magnitude of the average C-index in the validation cohort. The TCGA-LIHC cohort was used as the training set, and prognostic models were constructed from the 101 combinations. The concordance index (C-index) of each model was calculated in two validation cohorts (LIRI-JP and GSE116174). The model with the highest C-index was selected. Risk scores were computed as the sum of the product of gene expression and corresponding regression coefficients. Patients were classified into high-risk and low-risk groups based on the median risk score. Kaplan-Meier survival curves and receiver operating characteristic (ROC) curves were used to evaluate model performance.

Construction of nomograms
Nomograms were constructed by integrating risk scores with clinical variables, including age, gender, tissue grade, T-stage, M-stage, and N-stage. Calibration curves and ROC analyses were used to evaluate the predictive accuracy of the nomograms.

Functional enrichment analysis
Differentially expressed genes between the high-risk and low-risk groups (|log2FC| > 1, FDR < 0.05) were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses to explore their biological functions.

Immune status analysis
To assess immune status differences between the high-risk and low-risk groups, we calculated immune cell and immune function scores using single-sample gene set enrichment analysis (ssGSEA). Immune cell infiltration was estimated using seven algorithms: XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT. Correlations between immune cell infiltration and risk scores were evaluated using Pearson correlation.

Cell communication analysis
Single-cell transcriptomic data were normalized using the “NormalizeData” function, and the top 2,000 variably expressed genes were identified using the “FindVariableFeatures” function. The “CellChat” package in R was used to analyze cellular communication between tumor and immune cells.

Mutation analysis
Mutation analysis was performed using the “maftools” package to generate waterfall plots for the high-risk and low-risk groups. Microsatellite instability (MSI) differences between the two groups were analyzed, and survival analysis was conducted for high-MSI (H-MSI) and low-MSI (L-MSI) patients. Copy number variation (CNV) data for risk genes were also analyzed to assess the frequency of gene amplifications and deletions.

Analysis of immunotherapy and chemotherapy
Tumor immune dysfunction and exclusion (TIDE) scores were calculated for the high-risk and low-risk groups to predict their response to immunotherapy. The IMvigor210 cohort was used to evaluate the impact of anti-PD-L1 therapy on these groups. Chemotherapy sensitivity to four drugs (cisplatin, doxorubicin, epothilone B, and mitomycin C) was assessed using the “pRRophetic” package. Small molecule compounds with potential therapeutic effects were identified from the CMAP database by analyzing differentially expressed genes between the risk groups.

Identification of key genes
To further screen prognostic PTM-related genes, we employed three machine learning algorithms: Random Survival Forest (RSF), Least Absolute Shrinkage and Selection Operator (LASSO), and Support Vector Machine-Recursive Feature Elimination (SVM-RFE). Genes identified concurrently by all three algorithms were defined as key genes.

Cell culture
The cell lines used in this study (HL7702, Huh-7, HepG2, Hep3B, and SMMC7721) were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin at 37 °C in a humidified incubator with 5% CO₂. Cells were passaged every 2–3 days, and those in the logarithmic growth phase were used for subsequent experiments.

Quantitative Real-Time PCR (qRT-PCR)
Total RNA was extracted using TRIzol reagent and reverse-transcribed into cDNA using a reverse transcription kit. qPCR was performed using SYBR Green Mix. Primer sequences are listed in Table 1. Relative mRNA expression levels were calculated using the 2^−ΔΔCt method, with GAPDH serving as the internal control.

Western blot analysis
Total protein was extracted from cells using lysis buffer and quantified using the BCA protein assay. Equal amounts of protein were separated by SDS-PAGE and transferred to PVDF membranes. Membranes were blocked with 5% non-fat milk for 1 h, then incubated overnight at 4 °C with the primary antibody (21761-1-AP, 1:5000). After washing with TBST, membranes were incubated with the secondary antibody (SA00001, 1:5000) at room temperature for 1 h. Protein bands were visualized using ECL detection reagent, and band intensities were quantified using ImageJ software.

Cell proliferation and colony formation assays
For the CCK-8 assay, 5 × 10³ cells per well were seeded in 96-well plates. After incubation for 24 to 96 h, 90 µL of fresh culture medium and 10 µL of CCK-8 reagent were added to each well and incubated for 2 h. Absorbance was measured at 450 nm using a microplate reader.
For the colony formation assay, 1000 cells were seeded per well in 6-well plates and cultured for 10 days. Colonies were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. Colonies were counted using ImageJ software.

Transwell invasion and migration assays
For the invasion assay, the upper chamber of the Transwell insert was coated with Matrigel. A serum-free cell suspension (1 × 10⁵ cells) was added to the upper chamber, while the lower chamber contained medium supplemented with 20% FBS. After 48 h, non-invading cells were removed with a cotton swab. Invading cells were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and imaged under a microscope. Cell counts were analyzed using ImageJ. For the migration assay, the procedure was identical to the invasion assay except that the upper chamber was not coated with Matrigel.

Wound healing assay
Cells were seeded in 6-well plates and cultured until they reached approximately 90% confluency. A scratch was made using a pipette tip, and detached cells were removed by rinsing with PBS. The medium was then replaced with serum-free DMEM. The cells were incubated at 37 °C with 5% CO₂, and images of the wound area were captured at 0, 24, and 48 h using a microscope. The wound width was measured to evaluate cell migration.

Apoptosis assay
Cells were seeded in 6-well plates and cultured for 24 h. Apoptosis was assessed using an Annexin V-FITC Apoptosis Detection Kit, followed by flow cytometry analysis. Data were processed and analyzed using FlowJo software.

Statistical analysis
Bioinformatics analyses were performed using R software (version 4.3.2), with statistical comparisons between two groups conducted using the Wilcoxon test. Experimental data were analyzed and visualized using GraphPad Prism 10.3.1. Results are presented as mean ± standard deviation (SD). Comparisons between groups were assessed using Student’s t-test or one-way ANOVA. A p-value of < 0.05 was considered statistically significant.

Results

WGCNA network analysis
We collected genes related to seven types of PTMs—phosphorylation, acetylation, methylation, SUMOylation, Neddylation, ubiquitination, and glycosylation—from the MSigDB database. After removing duplicates, a total of 4,210 unique PTM-related genes were retained. WGCNA was performed on these genes, with a soft-thresholding power of 5 (Fig. 2A), resulting in the identification of 11 distinct modules (Figs. 2B–C). According to the module-trait relationship heatmap, the blue module showed the strongest correlation with cancer (Fig. 2D). Differential expression analysis of PTM-related genes, using a threshold of |log2FC| > 1.5 and FDR < 0.05, identified 533 DEGs, including 22 downregulated and 511 upregulated genes (Fig. 2E). Of these, 98 DEGs belonged to the blue module (Fig. 2F). Univariate Cox regression analysis of the 98 DEGs, with a significance threshold of p < 0.001, identified 23 prognostic genes (Fig. 2G).

Clustering analysis
Unsupervised clustering of HCC samples based on the expression of the 23 prognostic genes identified K = 2 determined as the optimal number of clusters (Fig. 3A). Heatmap analysis revealed that these prognostic genes were highly expressed in Cluster A and lowly expressed in Cluster B (Fig. 3B). Kaplan–Meier survival analysis showed that patients in Cluster A had significantly worse overall survival (OS) than those in Cluster B (Fig. 3C). These findings suggest that high PTM gene expression is associated with poor prognosis. Principal Component Analysis (PCA) further confirmed a clear distinction between the two clusters (Fig. 3D). Gene Set Variation Analysis (GSVA) revealed that Cluster A and Cluster B were significantly enriched in distinct biological pathways, including DNA replication, cell cycle, and homologous recombination (Fig. 3E). To investigate differences in immune-related functions, we also performed GSVA on immune-related pathways. The results indicated significant differences in T cell receptor signaling, Toll-like receptor signaling, and the complement and coagulation cascades between the two clusters (Fig. 3F).

Construction of the prognostic model
A prognostic model was constructed using the 23 identified prognostic genes (Table 2). The TCGA-LIHC cohort was used as the training set. A total of 101 combinations of machine learning algorithms were evaluated to develop the optimal prognostic model. The combination of Stepwise Cox regression (both directions) and Supervised Principal Components (StepCox[both] + SuperPC) yielded the highest average concordance index (C-index = 0.661) across two validation cohorts and was selected as the final model (Fig. 4A). Kaplan–Meier survival analysis showed that patients in the high-risk group had significantly lower OS than those in the low-risk group in the training and validation cohorts (Fig. 4B, D and F). The area under the curve (AUC) values for predicting 1-, 2-, 3-, and 4-year OS in the TCGA-LIHC cohort were 0.777, 0.692, 0.665, and 0.680, respectively (Fig. 4C). In the LIRI-JP cohort, AUC values were 0.748, 0.678, 0.663, and 0.722 (Fig. 4E), while in the GSE116174 cohort, they were 0.678, 0.630, 0.678, and 0.646 (Fig. 4G). These results collectively demonstrate that the model provides reliable predictive outcomes for HCC patients.

Clinical feature analysis
To explore the association between risk scores and clinical features, we analyzed the TCGA-LIHC cohort. Significant differences were observed in histological grade, clinical stage, and T stage between the high- and low-risk groups (Figure S1A). Further comparisons showed that higher-grade tumors, advanced clinical stages, and higher T stages were associated with higher risk scores (Figures S1B–D), while age showed no significant correlation (Figure S1E). Univariate Cox regression analysis identified clinical stage, T stage, M stage, and risk score as significant prognostic factors (Figure S1F). Multivariate Cox regression analysis confirmed that the risk score was an independent prognostic factor (Figure S1G).

Construction of the nomogram
To enhance the predictive accuracy of the prognostic model, we constructed a nomogram incorporating age, sex, histological grade, T stage, M stage, N stage, and risk score (Fig. 5A). Calibration plots demonstrated that the nomogram accurately predicted 1-, 2-, and 3-year OS (Fig. 5B). Multivariate ROC curve analysis showed that the nomogram outperformed the individual risk score and other clinical features in predictive accuracy (Figs. 5C–E).

Functional enrichment analysis
To investigate the biological differences between the high-risk and low-risk groups, we performed KEGG and GO enrichment analyses on the DEGs. The KEGG analysis revealed significant differences in pathways such as the cell cycle, Th1 and Th2 cell differentiation, and proteoglycans in cancer (Fig. 6A). GO analysis showed enrichment in pathways related to T cell activation regulation, immunoglobulin complex formation, and MHC class II receptor activity (Fig. 6B).

Immune status analysis
We applied ssGSEA to quantify immune cell infiltration and immune function in both risk groups. The results showed that scores for activated dendritic cells (aDCs), follicular helper T cells (Tfh), and regulatory T cells (Tregs) were higher in the high-risk group. Similarly, immune functions such as antigen-presenting cell co-inhibition, chemokine activity, and human leukocyte antigen (HLA) expression were more active in the high-risk group (Fig. 7A). To further assess immune infiltration, we used seven algorithms (XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT) and calculated Pearson correlation coefficients between immune cell infiltration and risk scores. Most immune cells showed a positive correlation with risk scores (Fig. 7B), suggesting enhanced immune activity in the high-risk group. Immune subtyping revealed a higher proportion of the C4 subtype (lymphocyte depletion) and a lower proportion of the C1 subtype (wound healing) in the high-risk group (Fig. 7C). Furthermore, the expression levels of immune checkpoint genes were significantly elevated in the high-risk group (Fig. 7D).

Cell communication analysis
Given the heightened immune activity observed in the high-risk group, we further analyzed cell–cell communication between tumor cells and immune cells using single-cell transcriptomic data. First, we annotated cell types based on sample information (Fig. 8A), and then examined the expression of the nine risk genes at the single-cell level, finding that all were expressed in malignant cells (Fig. 8B). Using the “CellChat” package, we analyzed intercellular communication between tumor cells and immune cells. The results revealed that tumor cells in the high-risk group exhibited more active signaling interactions via the ADIPONECTIN, LIGHT, LEP, and FGF pathways (Figs. 8C–F).

Genomic mutation analysis
We used the “maftools” package to generate waterfall plots illustrating genomic mutations in both risk groups. The high-risk group exhibited a higher overall mutation burden. The top three most frequently mutated genes in the high-risk group were TP53, CTNNB1, and TTN (Fig. 9A), whereas in the low-risk group, the top mutated genes were CTNNB1, TTN, and MUC16 (Fig. 9B). MSI is a crucial biomarker in cancer immunotherapy, was also analyzed [27]. The high-risk group demonstrated significantly higher MSI scores (Fig. 9C). Based on MSI status, HCC patients were classified into high-MSI (H-MSI) and low-MSI (L-MSI) groups. Survival analysis indicated better OS in the L-MSI group (Fig. 9D). Moreover, when combining risk scores with MSI status, the group with low risk and L-MSI showed the best prognosis (Fig. 9E). CNV analysis revealed the deletion and amplification rates of the nine risk genes (Fig. 9F). Finally, we analyzed the correlation between the expression of nine model genes and CNVs. The results indicated that ANAPC7, B4GALT3, SAMD1, DNAJC6, and CTSA showed significant correlations with CNVs, with ANAPC7 exhibiting the strongest correlation (Figure S2).

Immunotherapy analysis
Previous studies have shown that patients with TP53 mutations respond better to immunotherapy [28]. In our study, the high-risk group had a higher TP53 mutation frequency. Additionally, immune checkpoint genes were more highly expressed in the high-risk group. These findings suggest that high-risk patients may benefit more from immunotherapy. To evaluate immunotherapy response potential, we first calculated TIDE scores for both risk groups. The high-risk group had lower TIDE scores, indicating a potentially better response to immunotherapy (Fig. 10A). Next, we analyzed anti–PD-L1 treatment responses in the IMvigor210 cohort. Patients who responded to treatment had higher risk scores (Fig. 10B), and those in the high-risk group showed better survival outcomes following treatment (Fig. 10C). The model achieved an AUC of 0.661 in predicting immunotherapy response (Fig. 10D), demonstrating moderate predictive accuracy.

Drug sensitivity analysis
We employed the “pRRophetic” package to estimate the half-maximal inhibitory concentration (IC50) of four chemotherapeutic agents—cisplatin, doxorubicin, epirubicin B, and mitomycin C—in both high- and low-risk groups. The results showed that the IC50 values for all four drugs were significantly lower in the high-risk group (Fig. 10E–H), indicating higher sensitivity to chemotherapy in this group. To identify potential therapeutic compounds, we used the CMAP database to analyze DEGs between the high- and low-risk groups. Based on a connectivity score threshold of less than − 80, 18 small-molecule compounds were identified as potential therapeutic candidates (Table 3).

Identification of key genes
To further refine prognostic PTM-related genes, we applied three machine learning algorithms—LASSO, RSF, and SVM-RFE—to the 23 PTM prognostic genes (Fig. 11A–C). Through this multi-algorithm approach, ANAPC7 and CTSA were identified as key risk genes (Fig. 11D).
ROC curve analysis was conducted to evaluate their diagnostic value. ANAPC7 showed an AUC of 0.986 (Fig. 11E), and CTSA had an AUC of 0.954 (Fig. 11F), suggesting that ANAPC7 may have greater diagnostic utility. Consequently, ANAPC7 was considered a key gene warranting further investigation in HCC.

Bioinformatics and experimental validation of ANAPC7
Bioinformatics analyses revealed that ANAPC7 is highly expressed in HCC tissues (Figure S3A), and its elevated expression is associated with poor prognosis (Figure S3B). Clinically, ANAPC7 expression was significantly higher in tumors of higher grade, advanced clinical stage, and elevated T stage (Figures S3C–E). Multivariate Cox regression analysis confirmed that ANAPC7 is an independent prognostic factor for HCC (Figure S3F). Experimental validation further supported these findings. Western blot analysis confirmed overexpression of ANAPC7 in HCC tissues compared to normal liver cells (Fig. 12A). CCK8 assays demonstrated that silencing ANAPC7 significantly suppressed proliferation in Huh-7 and Hep3B cells (Fig. 12B). Transwell assays showed reduced invasive and migratory capacities upon ANAPC7 knockdown (Figs. 12C–D), and apoptosis assays revealed increased apoptotic rates following ANAPC7 silencing (Fig. 12E).

Experimental validation of the SAMD1 gene
Given the high coefficient of SAMD1 in the prognostic model, we further explored its functional role in HCC. RT-qPCR analysis revealed high SAMD1 expression in HCC tissues (Fig. 13A). CCK8 assays showed that SAMD1 knockdown significantly inhibited the proliferation of HepG2 and Hep3B cells (Fig. 13B), and this finding was further supported by colony formation assays (Fig. 13C). Transwell invasion assays demonstrated reduced invasiveness following SAMD1 silencing (Fig. 13D), while wound healing assays indicated impaired migratory ability (Fig. 13E).

Molecular docking analysis
We employed AutoDock Vina to perform molecular docking of the 18 small-molecule compounds identified from the CMAP database with ANAPC7 (Table 4). Among these, Digoxin, Digitoxin, Proscillaridin, and XMD-885 exhibited the lowest binding energies of − 9.0, − 8.6, − 8.4, and − 7.5 kcal/mol, respectively (Figures S4A–D), suggesting their potential as therapeutic agents targeting ANAPC7.

Discussion
HCC is a major global health threat characterized by high heterogeneity, strong invasiveness, and poor prognosis [29]. Although treatment strategies such as surgical resection, chemotherapy, and immunotherapy have brought new hope to patients, the 5-year survival rate for HCC remains below 20% [30, 31]. Therefore, there is an urgent need to develop novel prognostic models to better predict patient outcomes and inform clinical decision-making. PTMs regulate protein function by modifying protein structures, and accumulating evidence suggests that aberrant PTM expression can contribute to cancer development and progression by disrupting protein activity [32]. Thus, PTM-related gene signatures offer promising avenues for early diagnosis and prognosis prediction in HCC. In this study, we constructed a prognostic model based on PTM-related genes to provide a new approach for predicting HCC prognosis.
A total of 4,210 PTM-related genes were obtained from the MSigDB database. Using WGCNA, we identified gene modules closely associated with HCC. Prognosis-related genes within these key modules were then selected using univariate Cox regression analysis. Subsequently, we constructed 101 machine learning model combinations and selected the StepCox[both] + SuperPC model based on the highest average C-index in the validation cohort. Kaplan–Meier and ROC curve analyses in the validation cohort demonstrated that the model possessed robust predictive power.
ANAPC7, also known as APC7, is a subunit of the anaphase-promoting complex (APC) and functions as an E3 ubiquitin ligase that plays a vital role in mitotic progression. Previous studies suggest that decreased APC7 expression may be associated with poor prognosis in breast cancer [33]. However, current knowledge on the role of APC7 in cancer remains limited, and its mechanisms in tumorigenesis are still unclear. B4GALT3 (β−1,4-galactosyltransferase 3), a member of the B4GALT family, is involved in poly-N-acetyllactosamine biosynthesis [34]. Overexpression of B4GALT3 has been shown to promote the proliferation, invasion, and migration of neuroblastoma cells, while its inhibition suppresses the proliferation and migration of bladder cancer cells [35, 36]. CTSA (Cathepsin A) is a lysosomal serine protease that helps degrade intracellular and extracellular waste to maintain cellular homeostasis. Experimental studies indicate that CTSA is overexpressed in HCC and associated with poor prognosis [37]. Another study demonstrated that CTSA suppresses the proliferation, invasion, and migration of prostate cancer cells by inhibiting the p38 MAPK pathway [38]. DNAJC6, a member of the DNAJ heat shock protein family, is involved in encoding auxin. It is overexpressed in HCC and linked to adverse clinical outcomes. Additionally, DNAJC6 can promote epithelial–mesenchymal transition (EMT) via the TGF-β signaling pathway, thereby facilitating HCC progression [39]. IKBKE (Inhibitor of Nuclear Factor Kappa-B Kinase Subunit Epsilon) is an atypical inflammatory kinase that drives the progression of various cancers. IKBKE promotes breast cancer invasion and metastasis by phosphorylating Snail, thus preventing its degradation by the E3 ligase β-TRCP1 [40]. Furthermore, the IKBKE inhibitor CYT387 has been shown to suppress glioblastoma progression by activating the Hippo pathway [41]. MCRS1 (Microspherule Protein 1) is involved in multiple biological processes, including cell proliferation, migration, and senescence, and is considered a potential oncogene [42]. MCRS1 overexpression promotes proliferation, migration, and invasion in gastric cancer and has also been shown to promote lung cancer growth via the miR-155–Rb1 axis [42, 43].SAMD1 (Sterile Alpha Motif Domain Containing 1) interacts with unmethylated CpG islands and regulates promoter activity of numerous genes. Recent studies show that SAMD1 is overexpressed in pancreatic cancer and correlates with poor prognosis. Knockdown of SAMD1 inhibits pancreatic cancer cell invasion and migration [44]. ST6GALNAC4 (ST6 N-Acetylgalactosaminide Alpha-2,6-Sialyltransferase 4), a member of the sialyltransferase family, catalyzes the transfer of sialic acid from CMP-sialic acid to galactose-containing substrates. Studies indicate that ST6GALNAC4 induces aberrant glycosylation, thereby promoting proliferation, invasion, and migration in HCC. It also contributes to HCC progression by recruiting tumor-associated macrophages [45].TRIM16 (Tripartite Motif Containing 16) is an E3 ubiquitin ligase that marks substrate proteins for degradation. Abnormal expression of TRIM16 has been implicated in the development of several cancers. TRIM16 enhances glycolysis through the NIK–SIX1 axis, thereby promoting pancreatic cancer metastasis [46]. Additionally, CircPTK2 has been shown to overcome cisplatin resistance in non-small cell lung cancer via the miR-942/TRIM16 axis [47]. The nine key PTM genes identified in this study (such as ANAPC7, B4GALT3, CTSA, etc.) play crucial roles in HCC prognosis and immune regulation. This finding aligns with the broader trend in current cancer biomarker research, where an increasing number of non-classical cancer genes (such as metabolic enzymes and PTM regulators) are being validated as key prognostic and immune-related biomarkers. For instance, a recent pan-cancer study demonstrated that ADSL, a metabolic enzyme whose activity may be PTM-regulated, serves as a novel biomarker for the immune microenvironment across multiple cancers, including HCC [48]. This indirectly supports our strategy of screening novel biomarkers and therapeutic targets from the perspective of PTM regulatory networks. Our work further confirms that PTM-associated gene signatures not only predict patient prognosis but also profoundly correlate with tumor immune characteristics, offering novel molecular insights into the immune evasion mechanisms of HCC.
We performed KEGG and GO enrichment analyses to explore the biological functional differences between the high-risk and low-risk groups. The KEGG analysis revealed significant differences in pathways such as the cell cycle, Th1 and Th2 cell differentiation, and proteoglycans in cancer. Similarly, the GO enrichment analysis indicated notable differences in pathways related to T cell activation regulation, immunoglobulin complexes, and MHC class II receptor activity. These findings suggest that the immune system plays a critical role in distinguishing between the high- and low-risk groups. To further investigate immune-related features, we used ssGSEA to calculate immune cell infiltration and immune function scores for both groups. The high-risk group exhibited higher scores in immune cell populations such as aDCs, Tfh cells, and Tregs, as well as in immune functions including antigen presentation co-inhibition, HLA expression, and chemokine signaling. Dendritic cells (DCs) are among the most effective antigen-presenting cells and are pivotal in both innate and adaptive immune responses. DC-based vaccines have shown considerable promise in cancer immunotherapy [49]. Tfh cells, a subset of CD4 + T cells, provide essential help to B cells in response to various pathogens. Recent studies have demonstrated that Tfh cells enhance anti-tumor immunity and improve the efficacy of immunotherapy, correlating with favorable clinical outcomes [50]. In contrast, Tregs are immunosuppressive cells that infiltrate the tumor microenvironment, suppress anti-tumor immune responses, and facilitate immune evasion. The presence of Tregs is strongly associated with tumor proliferation, invasion, and metastasis [51]. Furthermore, immune cell infiltration levels calculated using seven different methods were mostly positively correlated with the risk score, indicating a higher degree of immune activity in the high-risk group. Immune checkpoint gene analysis also revealed significantly higher expression levels of immune checkpoint genes in the high-risk group. In these analyses of immune infiltration, we observed an intriguing phenomenon: the high-risk group exhibited active immune states alongside a higher prevalence of immunosuppressive cells and elevated expression of immune checkpoint genes. Tumor immune phenotypes can be categorized into three types: immune-inflammatory, immune-exclusion, and immune-desert [52]. The immune-inflammatory type is also termed a “hot tumor,” while immune-exclusion and immune-desert types are collectively referred to as cold tumors. Patients in the high-risk group clearly exhibit hot tumors, whereas those in the low-risk group present cold tumors. Hot tumors, characterized by high immune activity and high immunosuppression, respond favorably to immunotherapy, whereas cold tumors exhibit the opposite response. Converting cold tumors into hot tumors to enable immunotherapy is a current research focus. Literature indicates that promoting T-cell activation through immunostimulants, oncolytic viruses, chemotherapy, and radiotherapy; increasing antigen-specific T-cells via cell therapies and cancer vaccines; and enhancing T-cell trafficking and infiltration using oncogenic pathway inhibitors and CXCR4 inhibitors can facilitate this transformation, thereby making cold tumors more responsive to immunotherapy [53].Given the elevated immune activity in the high-risk group, we further investigated intercellular communication between tumor cells and immune cells in this group. The results showed that tumor cells in the high-risk group exhibited more active interactions with immune cells via signaling pathways such as ADIPONECTIN, LIGHT, LEP, and FGF. The ADIPONECTIN signaling pathway regulates local immune responses by modulating cytokine release from adipose tissue, thereby affecting tumor cell growth. It has been shown to suppress the proliferation and polarization of M1 macrophages, reduce B cell generation, and lower T cell reactivity, while promoting the proliferation and activation of M2 macrophages [54]. The LIGHT signaling pathway plays a key role in tumor immune responses and remodeling of the tumor microenvironment. It has been shown to enhance the function of effector cells, promote CD8 + T cell infiltration into tumors, and facilitate the formation of tertiary lymphoid structures [55]. The LEP pathway is involved in immunological and inflammatory processes in tumors and influences tumor progression by inhibiting T cell responses and promoting macrophage polarization toward the M2 phenotype. LEP signaling also enhances myeloid-derived suppressor cells (MDSCs) through IFN-γ production, which in turn modulates PD-L1 expression, leading to T cell apoptosis [56]. The FGF signaling pathway plays essential roles in cell proliferation, development, and differentiation, and its dysregulation is closely linked to the progression of various cancers. It has been particularly implicated in the interactions between pancreatic cancer and the tumor microenvironment [57]. Due to the elevated expression of immune checkpoint genes, the high-risk group may be more responsive to immunotherapy. We calculated the TIDE scores for both groups and found that the high-risk group had significantly lower scores, suggesting a better response to immunotherapy. To validate this, we analyzed the IMvigor210 cohort, a phase II clinical study evaluating the efficacy of the PD-L1 inhibitor atezolizumab in advanced urothelial carcinoma [58]. Consistent with our predictions, patients in the high-risk group exhibited better treatment responses and higher overall survival rates. Additionally, we assessed the sensitivity of both groups to four chemotherapeutic agents and found that the high-risk group was more sensitive to all four drugs. Patients in the high-risk group exhibit low survival rates, and clinical feature analysis indicates that higher risk scores correlate with more advanced clinical stages, suggesting greater tumor malignancy in this group. Tumors with higher malignancy exhibit accelerated proliferation rates. Since all four chemotherapeutic agents target cells in mitotic division, we hypothesize this may explain the heightened sensitivity of high-risk patients to these treatments. Based on the CMAP database, we also identified 18 small-molecule compounds with potential therapeutic effects. The clinical translation potential of this study warrants attention. A major bottleneck in current clinical treatment for hepatocellular carcinoma lies in the lack of reliable biomarkers capable of effectively guiding patient stratification, treatment selection, and clinical trial design. Our constructed PTM risk model, combined with potential therapeutic small molecules (e.g., digoxin) screened from the CMAP database and validated through molecular docking with the key target ANAPC7, provides crucial insights for developing novel combination therapy strategies and designing precision clinical trials. Literature indicates that tumor marker-based targeted drug discovery is pivotal for advancing personalized therapy [59]. Thus, our model serves not only as a prognostic prediction tool but also holds promise for future translation into systems supporting clinical decision-making. For instance, it could identify high-risk patient cohorts likely to benefit from immunotherapy or specific chemotherapies, or provide alternative drug options for treatment-resistant patients.
In conclusion, this study constructed a prognostic model for HCC based on PTMs, providing a novel approach for predicting patient outcomes. However, the study has certain limitations. Limited clinical sample stratification by etiological factors. The model was developed using data from public databases, and further validation with extensive clinical data is required to confirm its accuracy and clinical utility. Additionally, we conducted in vitro validation only for the ANAPC7 and SAMD1 genes, lacking in vivo validation and exploration of the functional roles of other model genes. Finally, since the transcriptomic data in the TCGA database originate from tissue samples, whether the model can enable non-invasive detection and application via liquid biopsy (such as ctDNA) will be a key direction for our future research.

Conclusion
In conclusion, we developed a robust prognostic model using PTM-related genes to predict survival outcomes and therapeutic responses in HCC patients. Our findings underscore the importance of PTMs in cancer progression and offer new avenues for improving HCC patient management through personalized treatment strategies.

Supplementary Information