Development and Validation of a Six-Gene Signature of Myeloid Antigen Presentation Dysfunction Based on Single-Cell and Multi-Cohort Transcriptomics for Predicting Prognosis and Recurrence of Hepatocellular Carcinoma

Abstract

Background: Myeloid antigen-presentation dysfunction is a central but incompletely translated immune state in hepatocellular carcinoma (HCC). This study aimed to develop a compact prognostic and recurrence model directly anchored in single-cell-defined myeloid antigen-presentation loss. Methods: In this prognostic prediction-model development and validation study, GSE149614 was used as the single-cell discovery cohort. Tumor and normal myeloid cells were extracted, antigen-presentation (AP) module scores were calculated, and AP-loss candidate genes were identified using pseudobulk aggregation. Candidate genes were projected to TCGA-LIHC, GSE14520, and GSE76427 bulk transcriptomic cohorts. Survival-oriented model refinement was performed in GSE14520 using overall survival (OS) and recurrence-free survival (RFS), followed by validation with Kaplan-Meier analysis, Cox regression, fixed-time AUC, clinicopathological comparison, pathway scoring, targeted cell-cell communication analysis, and nomogram construction. Results: A total of 13,784 myeloid cells were analyzed, including 8,209 tumor-infiltrating myeloid cells and 5,575 normal myeloid cells. Tumor-associated myeloid cells showed lower AP scores than normal myeloid cells (P = 0.04054). A compact six-gene signature consisting of SMOX, CSF1, AQP9, FLNB, COL7A1, and MXI1 was established. In GSE14520, the signature predicted poor OS (HR = 1.94, 95% CI: 1.41-2.65, P = 3.73e-05) and poor RFS (HR = 1.60, 95% CI: 1.23-2.08, P = 4.63e-04). In TCGA-LIHC, it also predicted inferior OS (HR = 1.49, 95% CI: 1.12-1.99, P = 5.83e-03). High-risk tumors showed enhanced glycolysis, hypoxia, epithelial-mesenchymal transition (EMT), angiogenesis, and IL6-JAK-STAT3 activity, with reduced antigen-presentation and IFNG response signals. Conclusions: We report a compact single-cell-derived myeloid AP-loss signature for HCC prognosis and recurrence stratification. The model links clinical risk to a biologically interpretable myeloid dysfunction state and broader tumor microenvironment remodeling.

1. Introduction

Hepatocellular carcinoma (HCC) is a primary liver malignancy with persistently high morbidity and mortality worldwide and represents one of the leading causes of cancer-related death.^1,2 Despite the gradual improvement of early screening systems for liver cancer, continuous optimization of local treatments such as surgical resection and local ablation, and ongoing advancement of systemic therapies including targeted agents and combined immunotherapy, the overall long-term survival benefit of HCC patients remains unsatisfactory. High postoperative recurrence rate and pronounced tumor heterogeneity have become core challenges limiting clinical efficacy. ³ Existing clinical indicators are inadequate for risk stratification of prognosis and recurrence, failing to meet the requirements of precision diagnosis and treatment. Therefore, molecular biomarkers with high specificity and clinical applicability are urgently needed. ⁴

The tumor microenvironment (TME) plays a critical regulatory role in the occurrence, progression, invasion, metastasis, and immune escape of HCC. ⁵ As the predominant immune cell population in the TME, myeloid cells serve as a key hub linking inflammation, tissue remodeling, and antitumor immunity. Myeloid cells not only initiate specific antitumor immune responses through antigen processing and presentation but also readily undergo functional polarization under the induction of the TME, acquiring an immunosuppressive phenotype that contributes substantially to tumor immune escape.^6-8 Current studies on myeloid cells in HCC have mostly focused on the analysis of infiltration abundance while neglecting the core importance of cellular functional status. Inactivation of antigen presentation in tumor-associated myeloid cells is a key mechanism underlying the failure of effective antitumor immune responses.^7,8 Compared with simple abundance evaluation, research centered on this functional state can more accurately reveal the core regulatory principles of tumor immunity, and molecular signatures constructed from this state may have stronger mechanistic interpretability and clinical translational value.^9,10

The development of single-cell RNA sequencing (scRNA-seq) has provided a powerful tool for high-resolution dissection of functional heterogeneity among cell populations in the TME. It enables precise identification of dysfunctional states such as impaired antigen presentation in myeloid cells and their associated molecular signatures, offering a novel perspective for mining tumor molecular biomarkers. However, in current cancer informatics research, a significant bottleneck remains in translating single-cell-level functional discoveries into routine clinical applications. Most studies fail to effectively project cell-level molecular signatures onto bulk transcriptomic data. ⁹ In contrast, previous HCC prognostic gene models based on bulk data commonly suffer from issues including excessive gene numbers, purely statistically driven construction logic, limited cross-cohort stability, and lack of direct association with specific functional states of the TME.^11-13 These limitations lead to ambiguous biological significance and restricted clinical translational potential of such models.

Against this background, the present study takes the inactivation of myeloid cell antigen presentation in HCC as the core entry point. Based on a scRNA-seq discovery cohort, we analyzed the differences in antigen presentation function between tumor and normal myeloid cells and screened candidate genes closely related to this inactivated state. Using a pseudobulk aggregation strategy, we projected single-cell gene signatures onto bulk transcriptomic data from multiple cohorts including TCGA-LIHC, GSE14520, and GSE76427. By integrating survival and recurrence outcomes, we refined and validated the model, ultimately establishing a compact and stable 6-gene prognostic and recurrence model. We further systematically evaluated the clinical value and underlying molecular mechanisms of this model from multiple dimensions, including survival prognosis, recurrence risk, correlation with clinical characteristics, pathway activity regulation, targeted cell-cell communication, incremental predictive value, and nomogram construction. This study aims to provide a novel molecular tool for prognosis evaluation and postoperative recurrence risk stratification in HCC patients, as well as clues for subsequent mechanistic research and translational exploration targeting myeloid functional reprogramming, thereby offering theoretical support for individualized management and precision diagnosis and treatment of HCC.

2. Methods

2.1. Data Sources and Study Design

The GSE149614 single-cell RNA sequencing dataset was used as the discovery cohort. ¹⁴ GSE14520 was used for survival-oriented model refinement and internal evaluation, ¹⁵ TCGA-LIHC was used for external overall survival validation, ¹⁶ and GSE76427 was used as an additional bulk transcriptomic cohort for cross-cohort gene projection. ¹⁷ This prognostic prediction-model development and validation study was reported in accordance with the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD) guideline, ¹⁸ and the completed TRIPOD checklist is provided as Supplementary Table S1. As a secondary analysis of public datasets, no new participant accrual was performed; dataset-level dates and follow-up information were obtained from the original source publications or database records. The overall study design consisted of four steps: first, identifying the inactivated state of antigen presentation in tumor-associated myeloid cells at the single-cell level; second, screening candidate genes stably associated with this state using a pseudobulk strategy ¹⁹ ; third, projecting candidate genes onto bulk cohorts and performing model refinement in combination with prognostic information ²⁰ ; and finally, evaluating the clinical and biological significance of the model via survival analysis, clinical correlation, pathway disparity, targeted cell-cell communication analysis, and nomogram construction. The overall workflow is shown in Figure 1.

Figure 1.

Overall study design and analytical workflow. This figure summarizes the single-cell discovery cohort, bulk validation cohorts, AP-loss candidate screening, survival-oriented refinement, validation, and biological interpretation steps

2.2. Preprocessing of Bulk Cohorts and Expression Matrix Processing

For TCGA-LIHC, the expression matrix and clinical data were matched at the patient level after data retrieval. Matrix matching and cleaning were performed using R version 4.5.2 and the data.table package. After excluding samples with missing key clinical outcome data or unmatched records, a total of 371 tumor samples were retained for subsequent analyses. For GSE14520, probe-to-gene symbol mapping was performed using the GPL571 platform annotation file, and for GSE76427, gene symbol annotation was completed using the GPL10558 platform. When multiple probes corresponded to the same gene, the probe with the highest average expression across samples was retained using the principle implemented in collapseRows ²¹ to reduce redundancy caused by multiple probes mapping to a single gene. After filtering and standardization, 11,248 common genes were identified across the three bulk cohorts, which were used for cross-cohort projection and model construction.

2.3. Extraction of Myeloid Cells From Single-Cell Data and Antigen Presentation Module Scoring

After processing the raw metadata of GSE149614, a total of 71,915 cells from 21 samples and 10 patients were included. Cells annotated as Myeloid in cell type and derived from either Tumor or Normal tissues were extracted, yielding a final set of 13,784 myeloid cells, including 8,209 tumor-infiltrating myeloid cells and 5,575 normal myeloid cells. Expression data were extracted from the raw count matrix according to the target cell barcodes using data.table, and normalization was performed using the log1p(CPM) method. ²² Specifically, raw expression values were converted to counts per million using the total counts per cell as library size, followed by log1p transformation to mitigate the influence of extreme values.

The antigen presentation (AP) module comprised the MHC-II module (HLA-DRA, HLA-DRB1, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQB1, CD74) and the MHC-I/antigen processing module (B2M, HLA-A, HLA-B, HLA-C, TAP1, TAP2, PSMB8, PSMB9, NLRC5).^23,24 The AP score for each cell was defined as the column mean of normalized expression values of the above genes. ²⁵

2.4. Pseudobulk Construction and Candidate Gene Screening From Single-Cell Data

Pseudobulk aggregation of myeloid cells was performed on a per-sample basis, ¹⁹ and a total of 18 myeloid pseudobulk samples were constructed, including 10 tumor samples and 8 normal samples. In tumor myeloid pseudobulk samples, samples were divided into the AP-high group and the AP-low group using the median total AP score of 7.069 as the cutoff. Three types of indicators were calculated for each gene separately: first, the difference in normalized log-transformed mean expression between tumor myeloid and normal myeloid cells; second, the difference in normalized log-transformed mean expression between the AP-low group and the AP-high group; third, the Spearman correlation between gene expression and tumor myeloid AP score. Considering the limited number of myeloid pseudobulk samples, we used a pre-specified direction- and effect-size-aware relaxation strategy rather than relying only on strict multiplicity-adjusted thresholds. Candidate genes were initially required to be upregulated in AP-low myeloid pseudobulk samples and tumor myeloid samples and negatively correlated with tumor myeloid AP score. If fewer than 8 candidate genes were obtained under the strict criteria, relaxed standards were adopted, requiring a normalized log mean expression difference of at least 0.15 between AP-low and AP-high groups, a mean difference of at least 0.15 between tumor and normal myeloid samples, a negative correlation with AP score no higher than -0.35, and a correlation P value less than 0.10. A total of 190 AP-loss candidate genes were finally obtained. In further prioritized feature screening, only genes with consistent upregulation in the AP-low group, stronger negative correlation with AP score, and an upregulation trend in tumor myeloid cells were retained, resulting in 19 prioritized signature genes for subsequent bulk projection and model refinement.

2.5. Construction of the 6-Gene Compact Model

In tumor samples from GSE14520, LASSO-Cox models were constructed for overall survival and recurrence-free survival separately using the 19 prioritized signature genes as input. ²⁶ Model fitting used the glmnet package with family set to cox, alpha set to 1,10-fold cross-validation, and a fixed random seed of 123.^27,28 The exact cross-validated tuning parameters were lambda.min = 0.09203 and lambda.1se = 0.17651 for OS and lambda.min = 0.11321 and lambda.1se = 0.16425 for RFS. Because the lambda.1se models retained no nonzero genes and lambda.min retained only one nonzero gene for each endpoint, we followed the pre-specified fallback strategy and selected the top 6 genes by the mean absolute coefficient values from the OS and RFS models, followed by recalibration according to single-gene Cox prognostic direction. Finally, 6 genes, SMOX, CSF1, AQP9, FLNB, COL7A1, and MXI1, were retained to form the compact risk model. The risk score was calculated as the weighted sum of standardized gene expression divided by the sum of absolute weights. The final formula was (0.4919 x SMOX - 0.1786 x CSF1 - 0.1680 x AQP9 + 0.1548 x FLNB - 0.0715 x COL7A1 + 0.0327 x MXI1)/1.0973. The refinement process from the 19-gene prioritized set to the 6-gene model is shown in Figure 3, and the final model weights are presented in Figure 4 and Supplementary Table S2.

2.6. Survival Analysis, Multivariate Analysis and Fixed Time-point AUC

In each cohort, the median value of the 6-gene risk score was used to define high-risk and low-risk groups. Kaplan-Meier curves were applied to evaluate overall survival and recurrence-free survival in GSE14520, as well as differences in overall survival in TCGA-LIHC. The survival and multivariable models were implemented with the survival package. The multivariate Cox model for GSE14520 included age, gender, tumor size, multiple nodules, AFP level and TNM stage; the multivariate model for TCGA-LIHC included age, gender, histological grade and clinical stage. Fixed time-point AUC was estimated using an event-based classification approach: samples with events occurring before a given time point were defined as the event group, the remaining samples as the non-event group, and samples censored before that time point were excluded. Overall survival in GSE14520 was evaluated at 12, 36 and 60 months, and recurrence-free survival at 12, 24 and 36 months; overall survival in TCGA-LIHC was evaluated at 1, 3 and 5 years.

2.7. Clinical Comparison, Pathway Scoring and Targeted Cell-Cell Communication Analysis

In the comparison of clinical characteristics between high- and low-risk groups, the Wilcoxon rank-sum test was used for continuous variables, ²⁹ and the chi-square test was used for categorical variables. Pathway analysis was performed using manually curated functional gene sets, and module scores were calculated for Hypoxia, Glycolysis, EMT, Inflammation, IL6_JAK_STAT3, Angiogenesis, IFNG_Response, and Antigen_Presentation. The module score was defined as the column mean of row-standardized expression values of genes within the module. Between-group comparisons were performed using the Wilcoxon test, and multiple comparisons were adjusted by the Benjamini-Hochberg method. Cell-cell communication analysis employed a targeted ligand-receptor scoring strategy instead of a panoramic framework such as CellChat or CellPhoneDB.^30,31 The sender populations were restricted to Tumor|Myeloid and Normal|Myeloid, and the receiver populations were restricted to Myeloid, T/NK, and Hepatocyte populations, with a focus on the CSF1-CSF1R, SPP1-CD44, SPP1-ITGAV/ITGB1, CCL2-CCR2, MIF-CD74, MIF-CXCR4, TNF-TNFRSF1A, TNF-TNFRSF1B axes, as well as HLA-related antigen presentation axes. The interaction score was defined as the square root of the product of the mean expression of the ligand and the mean expression of the receptor, ³² and the support score was defined as the square root of the product of the proportion of ligand-expressing cells and the proportion of receptor-expressing cells.

2.8. Comparison With Classic Clinical Indicators and Evaluation of Incremental Predictive Value

To evaluate the additional value of the 6-gene model compared with classic prognostic assessment methods for HCC, separate models were constructed in the GSE14520 cohort: AFP-only model, TNM-only model, AFP+TNM combined model, 6-gene-only model, and 6-gene+AFP+TNM combined model. In the TCGA-LIHC cohort, models were constructed as follows: Stage-only model, Stage+Grade model, 6-gene-only model, 6-gene+Stage model, and 6-gene+Stage+Grade model. Model performance was compared using the C-index and fixed time-point AUC. For GSE14520, overall survival was evaluated at 12, 36, and 60 months, and recurrence-free survival at 12, 24, and 36 months. For TCGA-LIHC, overall survival was assessed at 1, 3, and 5 years.

2.9. Construction of Combined Model With Classic Clinical Indicators, Nomogram and Calibration Curves

In the GSE14520 cohort, a combined Cox proportional hazards model was constructed based on the 6-gene risk score, AFP stratification and TNM stage, upon which a nomogram was established to predict 1-, 3- and 5-year overall survival probabilities. In the TCGA-LIHC cohort, the combined performance of the 6-gene score with clinical stage and histological grade was further evaluated. The nomogram model was built using the cph and nomogram functions in the rms package. Model calibration was performed using the bootstrap resampling strategy, and 1-, 3- and 5-year calibration curves were plotted to assess the consistency between predicted probabilities and actual observed outcomes.

2.10. Statistical Analysis

All data cleaning, statistical analysis, and plotting were performed in R version 4.5.2, with major analyses implemented using data.table, survival, glmnet, rms, ggplot2, patchwork, gridExtra, and scales. Correlation analyses were conducted using Spearman's rank correlation. This rank-based method was selected a priori for evaluating monotonic relationships between gene expression and module scores and for robustness in the small myeloid pseudobulk dataset; it was not selected after a formal normality test. Accordingly, no normality-test-based gatekeeping procedure was used to choose between Pearson and Spearman correlation. Survival differences were compared using the Kaplan-Meier method with the log-rank test, and multivariate analysis was performed using the Cox proportional hazards model. Unless otherwise specified, a two-sided P < 0.05 was considered statistically significant.

3. Results

3.1. Single-Cell Analysis Identifies Impaired Antigen Presentation Status in Tumor-Associated Myeloid Cells

In GSE149614, a total of 71,915 cells were collected from 21 samples and 10 patients. After extracting myeloid cells from tumor and normal tissues according to cell annotations, 13,784 myeloid cells were included for analysis, comprising 8,209 tumor myeloid cells and 5,575 normal myeloid cells. After constructing the AP score based on MHC-II and MHC-I/antigen processing-related genes, the AP score of tumor-associated myeloid cells was significantly lower than that of normal myeloid cells (P = 0.04054), suggesting impaired antigen presentation function in tumor myeloid cells (Figure 2).

Figure 2.

Identification of the myeloid antigen-presentation-loss state in hepatocellular carcinoma during single-cell discovery. Shown are differences in AP scores between myeloid cells derived from tumor and normal tissues in GSE149614, together with an overview of single-cell discovery of AP-loss candidate genes

3.2. Screening of AP-Loss Candidate Genes and Identification of the 19-Gene Prioritized Signature

At the single-cell level, a total of 10,845 genes were retained for subsequent analysis, and 18 myeloid pseudobulk samples were constructed, including 10 tumor samples and 8 normal samples. According to the preset combined screening criteria integrating the differences between AP-low and AP-high groups, the differences between tumor and normal myeloid cells, and the negative correlation with AP score, 190 AP-loss candidate genes were identified. In the prioritized feature screening with stricter emphasis on directional consistency and effect magnitude, 19 genes were finally retained as the core input for translating single-cell functional status into a clinical bulk model Figure 3 (Figures 2 and 3).

Figure 3.

Survival-oriented reduction workflow from the 19-gene candidate set to the compact 6-gene model. The figure shows the final compact genes and exact LASSO tuning parameters.

3.3. Survival-Driven Refinement Establishes the 6-Gene Compact Model

In the GSE14520 cohort, LASSO-Cox models for overall survival and recurrence-free survival were separately constructed based on the 19 prioritized signature genes, and the weights were recalibrated in accordance with the prognostic direction of individual genes. A 6-gene model consisting of SMOX, CSF1, AQP9, FLNB, COL7A1, and MXI1 was ultimately obtained. The exact OS model tuning parameters were lambda.min = 0.09203 and lambda.1se = 0.17651, and the RFS parameters were lambda.min = 0.11321 and lambda.1se = 0.16425. Final recalibrated weights were 0.4919 for SMOX, -0.1786 for CSF1, -0.1680 for AQP9, 0.1548 for FLNB, -0.0715 for COL7A1, and 0.0327 for MXI1 Figure 4. Compared with the initial 19-gene set, this model maintains the biological core derived from single-cell data while further improving compactness and interpretability. The refinement process from the 19-gene prioritized set to the 6-gene model is shown in Figure 3, and the final model weights are presented in Figure 4 and Supplementary Table S2.

Figure 4.

Composition and weights of the compact 6-gene model. Shown are the six genes included in the final model and their directions and relative weights in the risk score

To clarify the cellular context of the final genes, we further summarized their single-cell expression across major cell types and in myeloid pseudobulk samples. CSF1, SMOX, MXI1, AQP9, and FLNB were detectable in myeloid cells and showed variable expression across AP-low and AP-high tumor myeloid pseudobulk samples. At the same time, FLNB, AQP9, and COL7A1 also showed expression in non-myeloid compartments, including endothelial, fibroblast, and hepatocyte populations. These results indicate that the final signature should be interpreted as a myeloid AP-loss-derived risk program that captures both myeloid dysfunction and broader TME remodeling rather than as a purely myeloid-restricted marker set (Supplementary Figure S4 and Supplementary Tables S3-S4).

3.4. The 6-Gene Model Significantly Predicts Poor Survival and Recurrence in GSE14520

In GSE14520, the 6-gene model was significantly associated with overall survival (HR = 1.94, 95% CI: 1.41-2.65, P = 3.73e-05), and high-risk patients exhibited poorer long-term survival. In the same cohort, the model was also significantly associated with recurrence-free survival (HR = 1.60, 95% CI: 1.23-2.08, P = 4.63e-04), indicating that this signature can not only identify patients at high risk of death but also distinguish those with a stronger tendency toward postoperative recurrence. The risk stratification results in GSE14520 are shown in Figure 5, and the Kaplan-Meier curves are presented in Figure 7.

Figure 5.

Risk stratification results in the GSE14520 cohort. Shown are the distributions of risk scores, survival status, and the heatmap of gene expression

3.5. External Cohort Validation Supports the Cross-Platform Prognostic Value of the Model

In TCGA-LIHC, the 6-gene model was also significantly associated with inferior overall survival (HR = 1.49, 95% CI: 1.12-1.99, P = 5.83e-03), indicating that the model is reproducible across different platforms and population backgrounds. Although GSE76427 was mainly used for expression projection and consistency support rather than full survival endpoint validation in this study, overall results suggested that the myeloid antigen presentation inactivation signal reflected by this model exhibits favorable cross-cohort stability. Risk stratification in TCGA-LIHC is shown in Figure 6, and the corresponding overall survival curve is presented in Figure 7.

Figure 6.

Risk stratification results in the TCGA-LIHC cohort. Shown are the distributions of risk scores, survival status, and the heatmap of 6-gene expression in the external cohort

Figure 7.

Kaplan-Meier survival analyses of the 6-gene model. Shown are Kaplan-Meier curves for overall survival in GSE14520, recurrence-free survival in GSE14520, and overall survival in TCGA-LIHC, with log-rank P values annotated in the plots

3.6. Multivariate Analysis Indicates the Model Has Independent Prognostic Significance in GSE14520

In GSE14520, after incorporating the 6-gene risk score together with age, gender, tumor size, multiple nodules, AFP level and TNM stage into the multivariate Cox model, the risk score remained independently associated with overall survival (HR = 1.35, 95% CI: 1.05-1.73, P = 0.019). In recurrence-free survival analysis, the 6-gene risk score also retained independent statistical significance (HR = 1.27, 95% CI: 1.03-1.56, P = 0.025). Because GSE14520 was used for survival-oriented model refinement, these results should be interpreted as internal training-cohort validation. In contrast, the score showed a marginal association trend in the multivariate model of the external TCGA-LIHC cohort (HR = 1.16, P = 0.109), suggesting that differences in clinical composition, detection platforms and event structures between cohorts may affect the effect size of the model (Figure 8).

Figure 8.

Multivariable Cox regression forest plot for the 6-gene model and key clinical variables. Shown are the hazard ratios and 95% confidence intervals for the risk score and major clinical variables in multivariable Cox models

3.7. Fixed Time-point AUC Reveals Moderately Stable Predictive Performance of the Model

In GSE14520, the AUC values of the 6-gene model for 12-, 36- and 60-month overall survival were 0.730, 0.655 and 0.594, respectively, and for 12-, 24- and 36-month recurrence-free survival were 0.682, 0.633 and 0.621, respectively. In TCGA-LIHC, the AUC values for 1-, 3- and 5-year overall survival were 0.629, 0.616 and 0.633, respectively. Collectively, these results indicate that the model confers moderate and relatively stable time-dependent predictive ability (Figure 9).

Figure 9.

AUC results at different prespecified time points. This figure shows the AUC performance of the 6-gene model at different follow-up time points in GSE14520 and TCGA-LIHC

3.8. Systematic Differences in Clinical Characteristics and Pathway Profiles Between High- and Low-Risk Groups

In GSE14520, the high-risk group was more likely to exhibit elevated AFP levels (P = 9.92e-10) and showed a significant difference in age distribution (P = 0.009). In TCGA-LIHC, the high-risk group displayed a markedly different histological grade distribution (P = 1.74e-4), suggesting that the model is associated with a more aggressive clinical phenotype. These imbalances also indicate that part of the risk score signal overlaps with established clinical aggressiveness indicators; therefore, multivariable adjustment and combined-model evaluation were essential for assessing incremental value. The corresponding baseline clinical comparisons are shown in Tables 1 and 2. Pathway analysis revealed elevated module scores for Glycolysis, Hypoxia, EMT, Angiogenesis and IL6_JAK_STAT3 in the high-risk group, along with downregulated Antigen_Presentation and IFNG_Response. These findings imply that the high-risk state is accompanied by enhanced metabolic reprogramming, aggravated malignant phenotype and deepened immunosuppression. The overall pathway differences are presented in Figure 10, the distribution of representative pathways in Figure 11, and the correlations between the risk score and inflammatory markers as well as malignant pathway modules in Figure S1.

Table 1.

Comparison of Baseline Clinical Characteristics Between the Low-Risk and High-Risk Groups in the GSE14520 Cohort

Variable	Low-risk group	High-risk group	P value
Age, median	52.0	49.0	0.009
Sex, n	Female: 14; Male: 107	Female: 17; Male: 104	0.700
Tumor diameter, >5 cm/≤5 cm	Large: 40; Small: 80	Large: 48; Small: 73	0.375
Multinodular tumors, n	No: 95; Yes: 26	No: 95; Yes: 26	1.000
Cirrhosis, n	No: 10; Yes: 111	No: 9; Yes: 112	1.000
TNM stage, n	I: 49; II: 44; III: 2; IIIA: 13; IIIB: 7; IIIC: 0	I: 47; II: 34; III: 1; IIIA: 16; IIIB: 8; IIIC: 4	0.313
BCLC stage, n	0: 10; A: 83; B: 11; C: 11	0: 10; A: 69; B: 13; C: 18	0.386
AFP, >300 ng/mL/≤300 ng/mL	High: 31; Low: 88	High: 79; Low: 40	9.92e-10

Note. Continuous variables are presented as medians; categorical variables are presented as the number of cases in each subgroup.

Table 2.

Comparison of Baseline Clinical Characteristics Between the Low-Risk and High-Risk Groups in the TCGA-LIHC Cohort

Variable	Low-risk group	High-risk group	P value
Age, median	64.0	59.5	0.016
Sex, n	Female: 52; Male: 134	Female: 69; Male: 116	0.071
Histologic grade, n	G1: 42; G2: 87; G3: 51; G4: 4	G1: 13; G2: 90; G3: 71; G4: 8	1.74e-04
AJCC stage, n	Stage I: 98; Stage II: 38; Stage III: 2; Stage IIIA: 31; Stage IIIB: 5; Stage IIIC: 1; Stage IV: 1; Stage IVA: 0; Stage IVB: 1	Stage I: 73; Stage II: 48; Stage III: 1; Stage IIIA: 34; Stage IIIB: 3; Stage IIIC: 8; Stage IV: 1; Stage IVA: 1; Stage IVB: 1	0.147

Note. Continuous variables are presented as medians; categorical variables are presented as the number of cases in each subgroup.

Figure 10.

Heatmap of pathway differences between the low-risk and high-risk groups. Shown are the overall differences in pathways related to metabolism, inflammation, immunity, and malignant phenotypes

Figure 11.

Boxplots of representative differential pathways. Shown are differences in score distributions for representative functional pathways between the low-risk and high-risk groups

3.9. Targeted Cell-Cell Communication Analysis Supports Aberrant Signaling Output of Tumor-Associated Myeloid Cells

In the targeted cell-cell communication analysis, prominent activities were observed for axes including SPP1-CD44, SPP1-ITGAV/ITGB1, CCL2-CCR2, CSF1-CSF1R, MIF-CD74 and MIF-CXCR4 in the signaling output from Tumor|Myeloid to Tumor|T/NK and Tumor|Hepatocyte. Although several HLA-related antigen presentation axes retained moderate expression in both normal and tumor myeloid cells, the aberrant output of axes related to inflammatory recruitment and adhesion remodeling was more pronounced in tumor-associated myeloid cells. Because this analysis used a targeted ligand-receptor scoring strategy, the results should be interpreted as supportive mechanistic evidence rather than a comprehensive reconstruction of all intercellular communication events (Figure S2 and Figure S3).

3.10. The 6-Gene Model Provides Incremental Predictive Value Beyond Classic Clinical Indicators

In the overall survival analysis of GSE14520, the 6-gene model achieved an AUC of 0.730 at 12 months, which was higher than 0.642 for the AFP-only model and 0.666 for the TNM-only model. After combining the 6-gene score with AFP and TNM, the AUC values at 12, 36, and 60 months increased to 0.750, 0.735, and 0.722, respectively, with a C-index of 0.681, surpassing 0.643 for the AFP+TNM model and 0.623 for the 6-gene-only model. These results indicate that the model offers clear incremental predictive value over traditional clinical indicators (Table 3, Figures 12 and 13).

Table 3.

Comparison of Prognostic Performance Among Classic Clinical Indicators, the 6-Gene Model, and Combined Models

Cohort	Endpoint	Model	C-index	Time-dependent AUC
GSE14520	OS	AFP	0.582	12/36/60 months: 0.642/0.587/0.594
GSE14520	OS	TNM	0.626	12/36/60 months: 0.666/0.662/0.682
GSE14520	OS	6-gene	0.623	12/36/60 months: 0.730/0.655/0.594
GSE14520	OS	6-gene+AFP+TNM	0.681	12/36/60 months: 0.750/0.735/0.722
GSE14520	RFS	AFP	0.548	12/24/36 months: 0.587/0.568/0.563
GSE14520	RFS	TNM	0.575	12/24/36 months: 0.593/0.611/0.600
GSE14520	RFS	6-gene	0.602	12/24/36 months: 0.682/0.633/0.621
GSE14520	RFS	6-gene+AFP+TNM	0.629	12/24/36 months: 0.698/0.678/0.672
TCGA-LIHC	OS	Stage	0.591	1/3/5 years: 0.629/0.631/0.619
TCGA-LIHC	OS	Stage+Grade	0.610	1/3/5 years: 0.630/0.655/0.661
TCGA-LIHC	OS	6-gene	0.604	1/3/5 years: 0.629/0.616/0.633
TCGA-LIHC	OS	6-gene+Stage	0.639	1/3/5 years: 0.679/0.680/0.686
TCGA-LIHC	OS	6-gene+Stage+Grade	0.641	1/3/5 years: 0.676/0.688/0.696

Note. OS comparisons in the GSE14520 cohort included AFP, TNM, and their combined model; RFS comparisons used the same AFP-, TNM-, and AFP+TNM-based strategies. The TCGA-LIHC cohort included Stage, Stage+Grade, and models combining these clinical factors with the 6-gene signature. AUC values were estimated at prespecified time points.

Figure 12.

Comparison of C-index values between classic clinical indicators and the 6-gene model. Shown are changes in C-index across different cohorts and endpoints for AFP, TNM, Stage, and models combining these indicators with the 6-gene signature

Figure 13.

Comparison of time-dependent AUC values among classic clinical indicators, the 6-gene model, and combined models. Shown are the AUC changes at prespecified time points for overall survival in GSE14520, recurrence-free survival in GSE14520, and overall survival in TCGA-LIHC

In the recurrence-free survival analysis of GSE14520, the 6-gene model yielded AUC values of 0.682, 0.633, and 0.621 at 12, 24, and 36 months, respectively, all of which were superior to the AFP-only or TNM-only models. When combined with AFP and TNM, the corresponding AUC values further improved to 0.698, 0.678, and 0.672, and the C-index increased from 0.602 to 0.629. This demonstrates that the model can not only supplement traditional indicators in identifying mortality risk but also improve postoperative recurrence risk stratification (Table 3, Figures 12 and 13).

In the overall survival analysis of TCGA-LIHC, the clinical stage-only model showed AUC values of 0.629, 0.631, and 0.619 at 1, 3, and 5 years, whereas the 6-gene+Stage combined model increased these to 0.679, 0.680, and 0.686. With the further inclusion of Grade, the 3-year and 5-year AUC values of the combined model reached 0.688 and 0.696, respectively, with a C-index of 0.641, outperforming the Stage-only, Stage+Grade, and 6-gene-only models overall. This further supports that the signature provides complementary value to the classic HCC prognostic evaluation system (Table 3, Figures 12 and 13).

3.11. Combined Nomogram Further Improves the Clinical Translational Format of the Model

Based on the overall survival cohort of GSE14520, this study further constructed a nomogram model by integrating the 6-gene score, AFP level, and TNM stage. In this combined model, both the 6-gene score and TNM stage retained substantial contributions, with an HR of 1.364 for the 6-gene score and 3.314 for TNM stage III. The overall model C-index was 0.681, supporting its potential as a candidate tool for quantitative clinical risk assessment (Table 4, Figure 14).

Table 4.

Cox Regression Coefficients of the Nomogram Model for Overall Survival in the GSE14520 Cohort

Variable	HR	95% CI	P value
compact_signature_z	1.364	1.073-1.735	0.011
AFP_binHigh	1.135	0.722-1.785	0.583
TNM_advIII	3.314	2.113-5.195	1.78e-07

Note. The nomogram model was developed in the GSE14520 overall survival cohort and incorporated the standardized 6-gene risk score, AFP stratification, and TNM stage.

Figure 14.

Nomogram of the combined overall survival model in the GSE14520 cohort. The nomogram was built using the 6-gene score, AFP, and TNM stage to predict 1-, 3-, and 5-year overall survival probabilities

Calibration curves showed that the nomogram achieved favorable consistency in predicting 1-, 3-, and 5-year overall survival probabilities, with predicted values closely matching observed outcomes. This suggests that the combined model is not only statistically superior to traditional clinical models but also holds potential for clinical translation and application (Figure 15).

Figure 15.

Calibration curves of the combined overall survival model in the GSE14520 cohort. Shown is the calibration performance of the nomogram for predicting 1-, 3-, and 5-year overall survival, used to assess the agreement between predicted probabilities and observed outcomes

4. Discussion

This study centered on the inactivation of myeloid antigen presentation, extending from single-cell discoveries to multi-cohort bulk validation, and ultimately constructed a compact model consisting of 6 genes. Compared with statistical screening purely from bulk data, the advantage of this pipeline lies in its focused research theme, as the model is directly linked to a well-defined functional state in the TME, thereby facilitating a biologically interpretable research structure. Similar single-cell-to-bulk translational strategies have recently been used to derive HCC prognostic models from myeloid marker genes, SPP1-positive macrophage programs, and other immune-related cell states.^12,13 However, the present study differs by using impaired myeloid antigen presentation as the specific discovery anchor and by jointly evaluating prognosis, recurrence, clinical increment, and targeted intercellular signaling. Single-cell sequencing provides high-resolution support for dissecting cellular functional heterogeneity in the TME; however, translating single-cell-level functional findings into clinically applicable models stable across bulk cohorts remains a core challenge in cancer informatics. ⁹ In this study, pseudobulk aggregation was adopted to project single-cell functional genes to the sample level, ³³ followed by survival-driven LASSO-Cox model refinement. This approach preserved the core biological theme of myeloid antigen presentation inactivation while reducing the dimensional mismatch between single-cell and bulk data.

Another key advantage of the model lies in its compactness and dual-dimensional predictive value. In contrast to common prognostic models with dozens or even scores of genes in the field, the 6-gene model substantially reduces the difficulty of external validation and clinical translation costs, making it more convenient for result presentation, multi-center validation, and subsequent development of detection kits. More importantly, this model not only effectively predicts poor overall survival in HCC patients but also accurately identifies patients at high risk of postoperative recurrence, conferring more direct clinical implications for postoperative management of HCC. The high postoperative recurrence rate of HCC represents a critical bottleneck restricting long-term patient survival, ³⁴ and existing clinical indicators such as TNM stage and AFP show limited ability in recurrence risk stratification. ³⁵ The 6-gene model can serve as an important supplement to traditional clinical indicators, providing precise molecular evidence for stratified management and individualized intervention in patients at high postoperative recurrence risk. Meanwhile, this study further verified the incremental predictive value of the model: when combined with classic clinical indicators, it significantly improved the predictive performance for overall survival and recurrence-free survival in HCC patients. The nomogram constructed on this basis further transformed molecular signatures into an individualized quantitative risk tool directly applicable to clinical practice, strengthening the clinical translational potential of the model.

Mechanistically, the high-risk group exhibited enhanced pathways related to glycolysis, hypoxia, EMT, IL6-JAK-STAT3 signaling, and angiogenesis, accompanied by downregulated antigen presentation and IFNG response,^36,37 suggesting that the model reflects not a single molecular abnormality but an integrated malignant microenvironment state characterized by immunosuppression, inflammatory remodeling, and metabolic reprogramming. These results reveal that myeloid antigen presentation inactivation is not an isolated immune defect, but forms a vicious cycle with tumor metabolic abnormalities ³⁶ and malignant progression. ³⁸ Targeted cell-cell communication analysis further showed that tumor-associated myeloid cells transmit aberrant regulatory signals to surrounding T/NK cells, hepatocytes, and tumor cells via SPP1-, CCL2-, CSF1-, and MIF-mediated axes, thereby reinforcing the TME state featuring superimposed immunosuppression and malignant progression. ³⁹ In addition, the core genes in the 6-gene model each exert key functions in TME regulation: CSF1 modulates myeloid cell recruitment and immunosuppressive polarization through the CSF1-CSF1R axis, SMOX participates in polyamine metabolism, oxidative stress, and metabolic homeostasis, while COL7A1 and FLNB are associated with extracellular matrix remodeling, cytoskeletal regulation, and EMT-related phenotypes.^40-43 The single-cell expression summary further supports a dual interpretation: several genes are detectable in myeloid cells, whereas others also reflect non-myeloid stromal or hepatocyte-associated remodeling. The coordinated dysregulation of these genes collectively points to a malignant microenvironment related to myeloid antigen presentation inactivation in HCC.

The therapeutic implication of this finding is that patients with a high 6-gene score may represent a subgroup in which immune escape is coupled to myeloid inflammatory recruitment and matrix remodeling. Strategies aimed at myeloid reprogramming, including CSF1-CSF1R axis inhibition, TAM-targeted intervention, or rational combinations with immune checkpoint blockade, may therefore deserve further experimental evaluation in HCC models. ⁴⁰ The present model is not proposed as a treatment-selection biomarker at this stage, but it provides a candidate framework for linking transcriptomic risk stratification to future myeloid-directed therapeutic studies.

5. Limitations

This study has several limitations. First, all analyses were based on transcriptomic and clinical data from public databases, representing a retrospective informatics study lacking validation in an independent clinical cohort, and the association between core gene expression profiles and clinical phenotypes was not confirmed at the protein level. ⁴⁴ Second, GSE14520 was used both for survival-oriented refinement and for primary model evaluation; therefore, its performance should be interpreted as internal validation rather than a fully independent external validation. Third, heterogeneity exists among external validation cohorts in terms of detection platforms, ethnic backgrounds, and clinical follow-up information, which may weaken the consistent effect of the model, as reflected by the marginal association of independent prognostic value in the TCGA-LIHC cohort. Fourth, the cell-cell communication analysis employed a targeted ligand-receptor scoring strategy rather than a panoramic dedicated framework, rendering it more suitable as mechanistic support than conclusive evidence.^45,46 Future studies are warranted to validate the robustness of the model using independent clinical cohorts via qRT-PCR, immunohistochemistry, and other experimental methods. Meanwhile, in vitro and in vivo functional experiments are needed to clarify the specific molecular mechanisms by which core genes regulate myeloid antigen presentation. Further integration of spatial transcriptomics and single-cell functional assays will help refine the communication network between myeloid cells and other cell populations in the TME.^47,48

6. Conclusions

This study directly links HCC prognosis and recurrence stratification to the functional state of myeloid antigen presentation inactivation, differing from traditional research strategies focusing only on immune cell infiltration abundance. The constructed 6-gene model features clear biological interpretability and favorable clinical translational potential, serving as a novel molecular tool for prognosis evaluation and postoperative recurrence risk stratification in HCC patients, while providing important clues for subsequent mechanistic research and translational exploration targeting myeloid functional reprogramming. Overall, under the framework of single-cell functional state discovery, multi-cohort clinical projection, compact model validation, and clinical value expansion, this study offers new insights and strategies for transcriptomic informatics research and precision diagnosis and treatment of HCC.

Supplemental Material

Supplemental material - Development and Validation of a Six-Gene Signature of Myeloid Antigen Presentation Dysfunction Based on Single-Cell and Multi-Cohort Transcriptomics for Predicting Prognosis and Recurrence of Hepatocellular Carcinoma

Supplemental material for Development and Validation of a Six-Gene Signature of Myeloid Antigen Presentation Dysfunction Based on Single-Cell and Multi-Cohort Transcriptomics for Predicting Prognosis and Recurrence of Hepatocellular Carcinoma by Ye Tan, Jingxuan Xiang, Zehan Wang, Yu Wang, Aidong Chen, and Qiwen Wu in Cancer Informatics.

Supplemental Material

Supplemental material - Development and Validation of a Six-Gene Signature of Myeloid Antigen Presentation Dysfunction Based on Single-Cell and Multi-Cohort Transcriptomics for Predicting Prognosis and Recurrence of Hepatocellular Carcinoma

Supplemental material for Development and Validation of a Six-Gene Signature of Myeloid Antigen Presentation Dysfunction Based on Single-Cell and Multi-Cohort Transcriptomics for Predicting Prognosis and Recurrence of Hepatocellular Carcinoma by Ye Tan, Jingxuan Xiang, Zehan Wang, Yu Wang, Aidong Chen, and Qiwen Wu in Cancer Informatics.