Ferroptosis-mitochondria-related molecular subtypes and prognostic signature reveal immune and pathway heterogeneity in hepatocellular carcinoma

Abstract

BACKGROUND

Hepatocellular carcinoma (HCC) is highly heterogeneous, and conventional clinicopathological factors incompletely capture its biological diversity. Ferroptosis and mitochondrial dysfunction converge on oxidative stress, metabolism, and regulated cell death, but their integrated prognostic relevance remains unclear. We hypothesized that genes shared by ferroptosis and mitochondrial programs would define biologically distinct HCC subtypes and support prognostic stratification.

AIM

To characterize ferroptosis-mitochondria-related heterogeneity and develop a prognostic signature for HCC.

METHODS

Ferroptosis-related genes from FerrDb were intersected with mitochondrial genes from MitoCarta 3.0 and analyzed in The Cancer Genome Atlas liver HCC cohort. Differential expression, Cox regression, consensus clustering, immune scoring, pathway enrichment, least absolute shrinkage and selection operator Cox modeling, 7:3 internal validation, and external evaluation in GSE14520 and the International Cancer Genome Consortium Liver Cancer-RIKEN Japan (ICGC-LIRI-JP) cohort were performed. A nomogram was constructed and evaluated using calibration and decision curve analyses; exploratory analysis of public proteomic data was also performed.

RESULTS

Forty-seven overlapping genes were identified, significantly exceeding the overlap expected by chance (1.71-fold enrichment; P = 2.41 × 10^-4). Four molecular subtypes showed distinct survival and immune score patterns. A six-gene score comprising NT5DC2, TIMM9, CISD1, LONP1, PRDX6, and LIG3 was constructed. In the internal validation cohort, the model achieved 1-, 3-, and 5-year areas under the curve of 0.767, 0.740, and 0.575, respectively, with a C-index of 0.654. In the external ICGC-LIRI-JP cohort, the corresponding areas under the curve were 0.682, 0.648, and 0.581, respectively. In both GSE14520 and ICGC-LIRI-JP, the continuous risk score remained significantly associated with overall survival, whereas the dichotomized risk groups did not show statistically significant Kaplan-Meier separation. High-risk tumors showed proliferative and metabolic pathway enrichment.

CONCLUSION

Ferroptosis-mitochondria-related genes capture prognostic and immunological heterogeneity in HCC, although the six-gene score requires further prospective and experimental validation.

Key Words: Hepatocellular carcinoma; Ferroptosis; Mitochondria; Prognostic signature; Immune cell scores; Bioinformatics

Core Tip: This study integrated ferroptosis- and mitochondria-related genes to characterize hepatocellular carcinoma heterogeneity and develop a six-gene prognostic score. Four molecular subtypes showed distinct survival and immune score patterns. Internal validation supported short- and intermediate-term discrimination, whereas performance weakened at 5 years, and Kaplan-Meier separation after risk-group dichotomization was not statistically significant in the external cohorts. The principal contribution is a compact, biologically interpretable framework linking ferroptosis, mitochondrial dysfunction, immune remodeling, and metabolic pathway alterations in hepatocellular carcinoma.

INTRODUCTION

Hepatocellular carcinoma (HCC), the predominant histological subtype of primary liver cancer, remains a major global health burden. According to the most recent global cancer statistics, liver cancer continues to rank among the leading causes of cancer-related death worldwide[1]. Despite advances in surveillance, locoregional therapy, systemic treatment, and immunotherapy, clinical outcomes remain poor for many patients with HCC, largely because of late diagnosis, a high risk of recurrence, and marked molecular heterogeneity[2,3]. Conventional clinicopathological staging does not fully capture the biological heterogeneity underlying differences in prognosis and treatment response among patients at the same clinical stage. Therefore, biologically interpretable molecular classifications and prognostic biomarkers are still needed to improve risk stratification and facilitate translational research.

Ferroptosis is an iron-dependent form of regulated nonapoptotic cell death characterized by lethal lipid peroxidation[4,5]. It is closely linked to iron metabolism, redox imbalance, glutathione metabolism, and membrane lipid remodeling. Increasing evidence indicates that ferroptosis is involved in tumor development, therapy resistance, metabolic vulnerability, and stress adaptation[6,7]. Beyond its role in tumor cell death, ferroptosis may also influence immune surveillance and tumor microenvironment remodeling through lipid peroxidation, inflammatory signaling, and stress-related metabolic changes[6,7]. Because HCC arises in a metabolically active organ and is frequently accompanied by oxidative and inflammatory stress, ferroptosis-related mechanisms may be particularly relevant to HCC biology.

Mitochondria are central regulators of energy production, reactive oxygen species generation, metabolic reprogramming, apoptotic signaling, and cellular adaptation to stress. Mitochondrial dysfunction is a hallmark of malignant transformation and tumor progression[8,9]. Moreover, mitochondrial metabolism and redox regulation are closely linked to ferroptosis susceptibility. Mitochondria may modulate ferroptosis susceptibility by regulating reactive oxygen species production, iron-sulfur cluster homeostasis, lipid metabolism, and cellular energy stress[10,11]. Genes involved in both ferroptosis and mitochondrial biology may therefore reflect the intersection of oxidative cell death, metabolic stress adaptation, immune remodeling, and HCC progression.

Although ferroptosis-related signatures[12] and mitochondrial signatures[13] have been explored separately, most transcriptomic studies have treated these two gene programs as independent biomarker sources. This strategy may overlook genes at the interface between ferroptosis-associated oxidative stress and mitochondrial dysfunction. The integrated prognostic and immunological significance of ferroptosis-mitochondria-related genes in HCC therefore remains insufficiently defined. In this study, we intersected ferroptosis-related genes from FerrDb with mitochondrial genes from MitoCarta 3.0, identified ferroptosis-mitochondria-related genes, and systematically investigated their associations with molecular subtypes, prognosis, immune score patterns, pathway enrichment, and proteomic evidence from public datasets in HCC. Accordingly, we hypothesized that genes at the interface of ferroptosis and mitochondrial biology would define biologically distinct HCC subtypes and support prognostic stratification.

MATERIALS AND METHODS

Data acquisition and preprocessing

Transcriptomic and clinical data for HCC were obtained from The Cancer Genome Atlas liver HCC (TCGA-LIHC) cohort. External cohorts included GSE14520 from the Gene Expression Omnibus and LIRI-JP from the International Cancer Genome Consortium (ICGC). The Gene Expression Omnibus is a public repository for high-throughput functional genomics datasets[14], and ICGC provides large-scale cancer genomic resources[15,16]. The GSE14520 cohort has been widely used as an external HCC transcriptomic dataset with survival information[17]. Following data preprocessing, 366 patients from TCGA-LIHC, 221 from GSE14520, and 93 from ICGC-LIRI-JP were included in the survival analyses; all patients had matched expression and survival data. For TCGA-LIHC, gene expression data were analyzed as log2-transformed TPM values. For GSE14520, probe-level microarray data were mapped to gene symbols and summarized at the gene level. For ICGC-LIRI-JP, FPKM values were log2-transformed using log2(FPKM + 1).

Construction of the ferroptosis-mitochondria-related gene set

Ferroptosis-related genes were obtained from FerrDb, a manually curated database of ferroptosis regulators and ferroptosis-disease associations[18]. Mitochondrial genes were obtained from MitoCarta 3.0, which provides a curated inventory of human mitochondrial proteins with suborganelle localization and pathway annotations[19]. After duplicate gene symbols were removed, 484 ferroptosis-related genes and 1136 mitochondrial genes were included in the intersection analysis. The 47 overlapping genes were designated as ferroptosis-mitochondria-related genes and were included in subsequent differential expression, survival, molecular subtype, immune scoring, pathway enrichment, and prognostic modeling analyses.

To determine whether the overlap between the ferroptosis-related and mitochondrial gene sets was greater than expected by chance, a one-sided hypergeometric test was performed using 20000 genes as the reference background. The test incorporated 484 ferroptosis-related genes, 1136 mitochondrial genes, and the 47 observed overlapping genes. The expected overlap and fold enrichment were also calculated.

Differential expression and survival analysis

Differential expression analysis between HCC tissues and normal liver tissues in TCGA-LIHC was performed using the limma framework[20]. Genome-wide differential expression results were generated, after which the ferroptosis-mitochondria-related genes were examined separately. Among the 47 overlapping genes, significantly dysregulated genes were defined using an adjusted P value < 0.05 and an absolute log2 fold change > 1. Univariate Cox regression analysis was used to evaluate the association between each overlapping gene and overall survival. Genes meeting the nominal significance threshold of P < 0.05 and those remaining significant after false discovery rate (FDR) correction were recorded separately.

Consensus clustering and molecular subtype analysis

Consensus clustering was performed using the expression profiles of the 47 ferroptosis-mitochondria-related genes in TCGA-LIHC patients. The ConsensusClusterPlus algorithm was used to evaluate clustering stability for k values ranging from 2 to 6[21]. Resampling-based consensus clustering was conducted using a fixed random seed, with repeated subsampling of patients and features to assess cluster robustness. The molecular subtype solution was selected based on the consensus matrix, cumulative distribution function curve, delta area plot, cluster size, and survival differences. Kaplan-Meier analysis was used to compare overall survival among molecular subtypes.

Construction of the prognostic risk signature

Survival-associated ferroptosis-mitochondria-related genes were subjected to least absolute shrinkage and selection operator (LASSO) Cox regression analysis to construct a prognostic signature. LASSO regression was implemented using the glmnet framework[22] with a fixed random seed (20260520). Tenfold cross-validation was performed, and lambda.min was selected as the penalty parameter. The final risk score was calculated as the sum of each selected gene’s expression value multiplied by its corresponding Cox coefficient. Patients in TCGA-LIHC were stratified into high- and low-risk groups according to the median risk score. In the external cohorts, the risk score was analyzed both as a continuous variable in Cox regression and as a dichotomized variable for Kaplan-Meier visualization.

To address potential overfitting, an additional internal validation analysis was performed. The 366 TCGA-LIHC patients were randomly divided into a training cohort (n = 256) and an internal validation cohort (n = 110) at a 7:3 ratio using a fixed random seed (20260520). Expression values of the six model genes were standardized using the means and standard deviations derived exclusively from the training cohort. A multivariable Cox model containing the six predefined genes was fitted to the training cohort, and the resulting coefficients were used without refitting to calculate risk scores in the internal validation cohort. The median risk score in the training cohort was used as the cutoff for both cohorts. Kaplan-Meier analysis, univariate Cox regression using the continuous risk score, and time-dependent receiver operating characteristic (ROC) analysis at 1, 3, and 5 years were then performed separately in the training and internal validation cohorts. Harrell’s C-index was additionally calculated for each subset, with higher risk scores specified as indicating shorter survival, so that higher concordance values represented better discrimination.

External performance evaluation and independent prognostic analysis

The external performance of the risk score was assessed in the independent GSE14520 and ICGC-LIRI-JP cohorts. Kaplan-Meier survival analysis, univariate Cox regression, and time-dependent ROC analysis were performed. Time-dependent ROC curves were generated using the timeROC framework to evaluate 1-, 3-, and 5-year predictive performance[23]. After platform-specific preprocessing, expression values for each model gene were standardized within each external cohort using gene-wise z scores. The same TCGA-derived model gene set and coefficients were then used without model refitting. A cohort-specific median risk score was used only for dichotomized Kaplan-Meier visualization, whereas the continuous risk score was used for Cox and ROC analyses. Because the three cohorts were analyzed separately rather than pooled, ComBat adjustment was not applied; within-cohort standardization reduced scale differences but could not eliminate all residual platform or batch effects. Hazard ratios (HRs) used to summarize external performance were derived from univariate Cox regression models treating the risk score as a continuous variable. Multivariable Cox regression analysis was performed in TCGA-LIHC to determine whether the risk score was independently associated with overall survival after adjustment for age, sex, pathological stage, and tumor grade.

To evaluate the contribution of NT5DC2 to model performance, a leave-one-gene-out sensitivity analysis was performed by removing the NT5DC2 term from the six-gene score while retaining the original coefficients of TIMM9, CISD1, LONP1, PRDX6, and LIG3 without refitting. The resulting five-gene score was evaluated in TCGA-LIHC and GSE14520 using Kaplan-Meier analysis, univariate Cox regression using the continuous score, and 1-, 3-, and 5-year time-dependent ROC analysis.

Nomogram, calibration curve, and decision curve analysis

A prognostic nomogram was constructed by integrating the six-gene risk score with age, sex, pathological stage, and tumor grade. These clinicopathological variables were retained a priori for clinical adjustment and comparability, even when an individual coefficient was not statistically significant. A single multivariable Cox model was used to generate separate 1-, 3-, and 5-year overall survival probability scales. Nomogram construction and calibration were performed using the rms framework. Calibration curves were generated with bootstrap resampling under a fixed random seed (20260520; 200 resamples) to assess agreement between predicted and observed survival probabilities. Decision curve analysis was used to compare the potential clinical net benefit of the risk score-only, clinical, stage-only, and combined risk-plus-clinical models[24]. For holdout internal validation of discrimination, the same predefined 7:3 TCGA split was used; after exclusion of incomplete clinical records, the nomogram training and internal validation subsets contained 238 and 102 patients, respectively. The training-derived risk score and Cox coefficients were applied to the internal validation subset without refitting, and Harrell’s C-index was calculated for the risk score-only and combined models. After internal validation, the final nomogram displayed in the manuscript was refitted using all 340 TCGA patients with complete clinical information. Regression modeling procedures followed established survival modeling strategies[25].

Immune score analysis

Single-sample gene set enrichment analysis (ssGSEA) provides sample-level enrichment estimates and has been widely used to infer relative pathway or immune cell activity from bulk expression data[26]. In this study, ssGSEA was implemented using the gene set variation analysis (GSVA) framework to calculate relative immune cell scores for 25 immune cell populations based on previously reported immune cell gene signatures[27,28]. These scores represent transcriptome-derived enrichment estimates rather than direct cell counts or histological measurements. Differences among molecular subtypes and between risk groups were assessed using nonparametric tests. Spearman correlation analysis was used to evaluate associations between the continuous risk score and immune cell scores. Benjamini-Hochberg FDR correction was applied to the 25 immune cell comparisons, and FDR < 0.05 was considered statistically significant.

Functional pathway enrichment analysis

GSVA was used to estimate Hallmark pathway activity scores for each sample[29]. Gene set enrichment analysis (GSEA) was performed to identify pathways enriched between the high- and low-risk groups[30]. Hallmark gene sets were obtained from the Molecular Signatures Database[31]. Pathway analyses were based on genome-wide transcriptomic profiles and independent Hallmark gene sets rather than restricted to the six genes included in the prognostic signature. GSVA scores were compared between risk groups, and GSEA was performed using a preranked gene list derived from high- vs low-risk differential expression statistics. Pathways with adjusted P values below the significance threshold were considered significantly enriched.

Public CPTAC-based protein-level exploration

Public CPTAC proteomic data accessed through UALCAN were used to explore the protein abundance patterns of the six model genes in HCC[32]. The protein abundance of NT5DC2, TIMM9, CISD1, LONP1, PRDX6, and LIG3 was compared between normal liver tissues and primary HCC tissues. UALCAN-CPTAC results were used as exploratory protein-level observations from public data, not as in-house experimental validation. Additional protein-level information was obtained from the Human Protein Atlas[33].

Statistical analysis

All statistical analyses were performed using R software. The principal R packages used in this study included limma, ConsensusClusterPlus, glmnet, survival, survminer, timeROC, rms, GSVA, clusterProfiler[34], and standard visualization packages. Wilcoxon rank-sum tests or Kruskal-Wallis tests were used for group comparisons as appropriate. Survival differences were evaluated using Kaplan-Meier curves and log-rank tests. Cox proportional hazards models were fitted to evaluate associations with overall survival. P < 0.05 was considered statistically significant unless otherwise specified. FDR correction was applied to account for multiple testing where appropriate.

RESULTS

Identification of ferroptosis-mitochondria-related genes in HCC

To systematically investigate the role of ferroptosis-mitochondria-related genes in HCC, ferroptosis-related genes were obtained from FerrDb, and mitochondrial genes were obtained from MitoCarta 3.0. The intersection of the two gene sets yielded 47 ferroptosis-mitochondria-related genes, which were included in subsequent analyses (Figure 1A and B).

Figure 1 Identification of ferroptosis-mitochondria-related genes in hepatocellular carcinoma. A: Study workflow; B: Venn diagram showing the overlap between ferroptosis-related and mitochondrial genes; C: Expression heatmap of the 12 significantly differentially expressed ferroptosis-mitochondria-related genes in tumor and normal liver tissues from The Cancer Genome Atlas liver hepatocellular carcinoma cohort; D: Differential expression volcano plot; E: Forest plot of representative survival-associated overlapping genes based on univariate Cox regression. TCGA-LIHC: The Cancer Genome Atlas liver hepatocellular carcinoma.

The statistical significance of this overlap was further evaluated using a one-sided hypergeometric test. Based on a reference background of 20000 genes, the expected overlap between the two gene sets was 27.49 genes, whereas 47 overlapping genes were observed, corresponding to a 1.71-fold enrichment (P = 2.41 × 10^-4). These results indicate that the overlap between ferroptosis-related and mitochondrial genes was significantly greater than expected by chance.

We first performed genome-wide differential expression analysis between TCGA-LIHC tumor and normal tissues, and then focused on the 47 ferroptosis-mitochondria-related genes. Among these genes, 12 were significantly differentially expressed. In parallel, univariate Cox regression analysis was performed to evaluate the associations between the 47 genes and overall survival. Seventeen genes were nominally associated with overall survival at P < 0.05, of which ten remained significant after FDR correction. Integrated analysis identified NT5DC2, PDK4, and GSTZ1 as both differentially expressed and associated with overall survival; among them, NT5DC2 also met the FDR-adjusted survival significance criterion (Figure 1C-E).

Consensus clustering identified four ferroptosis-mitochondria-related molecular subtypes

To determine whether ferroptosis-mitochondria-related genes could define molecular heterogeneity in HCC, consensus clustering was conducted using the expression profiles of the 47 overlapping genes in TCGA-LIHC patients. Based on the consensus matrix and cluster stability assessment, k = 4 was selected as the optimal clustering solution. Patients were assigned to four molecular subtypes: C1 (n = 68), C2 (n = 135), C3 (n = 58), and C4 (n = 105) (Figure 2A-D).

Figure 2 Identification of ferroptosis-mitochondria-related molecular subtypes in The Cancer Genome Atlas liver hepatocellular carcinoma. A: Consensus clustering cumulative distribution function curve; B: Delta area plot; C: Consensus matrix for k = 4; D: Heatmap showing the expression patterns of the 47 genes across the four molecular subtypes; E: Kaplan-Meier survival curves comparing overall survival among the four molecular subtypes; F: Distribution of clinical characteristics across the four molecular subtypes. TCGA-LIHC: The Cancer Genome Atlas liver hepatocellular carcinoma; CDF: Cumulative distribution function.

Kaplan-Meier survival analysis showed significant differences in overall survival among the four molecular subtypes (log-rank P = 0.0006221), indicating that ferroptosis-mitochondria-related transcriptional patterns were associated with HCC prognosis (Figure 2E). The distributions of age, sex, pathological stage, and tumor grade were further visualized across the four subtypes, showing subtype-specific patterns in clinicopathological characteristics (Figure 2F). These findings suggest that ferroptosis-mitochondria-related genes may help stratify patients with HCC into biologically and clinically distinct subgroups.

Molecular subtypes displayed distinct ssGSEA-derived immune cell score patterns

To further characterize biological differences among the four molecular subtypes, ssGSEA was used to calculate relative scores for 25 immune cell populations. All 366 TCGA-LIHC patients with molecular subtype labels were included. Kruskal-Wallis tests followed by Benjamini-Hochberg correction demonstrated broad immune heterogeneity: 23 of 25 immune cell scores differed significantly among the four subtypes after multiple-testing correction (FDR < 0.05).

Notably, C3 showed the lowest median scores for most immune cell populations, including plasmacytoid dendritic cells, activated CD4 T cells, natural killer T cells, activated CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, natural killer cells, regulatory T cells, neutrophils, activated dendritic cells, activated B cells, eosinophils, monocytes, memory B cells, and plasma cells. In contrast, C1 showed relatively higher scores for several T-cell and natural killer cell-related populations, whereas C4 showed higher scores for immature B cells, dendritic cells, eosinophils, monocytes, and memory B cells. These results are consistent with an immune-low pattern in C3. However, because no established immune-hot/immune-cold classification framework was directly applied, C3 should not be considered equivalent to a formally defined immune-cold HCC class (Figure 3).

Figure 3 Single-sample gene set enrichment analysis-derived immune cell score patterns among the four molecular subtypes. A: Mean scores of 25 immune cell populations across the four molecular subtypes; B: Boxplots with individual patient observations for representative immune cell scores showing significant subtype differences after Benjamini-Hochberg correction; C: Sample-level heatmap of representative immune cell scores ordered by molecular subtype. ssGSEA: Single-sample gene set enrichment analysis; FDR: False discovery rate.

Construction of a six-gene ferroptosis-mitochondria-related prognostic signature

To establish a prognostic model based on ferroptosis-mitochondria-related genes, survival-associated candidate genes were subjected to LASSO Cox regression analysis in the TCGA-LIHC cohort. A six-gene prognostic signature comprising NT5DC2, TIMM9, CISD1, LONP1, PRDX6, and LIG3 was subsequently derived (Figure 4A and B). The risk score was calculated using the following formula: Risk score = 0.211 × NT5DC2 + 0.149 × TIMM9 + 0.150 × CISD1 + 0.145 × LONP1 + 0.201 × PRDX6 + 0.0485 × LIG3. Coefficients are rounded for presentation; full-precision coefficients were used in all analyses.

Figure 4 Construction and internal validation of the six-gene ferroptosis-mitochondria-related prognostic signature. A: Least absolute shrinkage and selection operator coefficient profiles; B: Tenfold cross-validation curve for selection of the optimal penalty parameter; C: Risk score distribution in the full The Cancer Genome Atlas liver hepatocellular carcinoma (TCGA-LIHC) cohort; D: Survival status distribution in the full cohort; E: Heatmap of the six model genes ordered by risk score; F: Kaplan-Meier overall survival curves for the full TCGA-LIHC high- and low-risk groups; G: Time-dependent receiver operating characteristic curves for 1-, 3-, and 5-year overall survival in the full TCGA-LIHC cohort; H: Kaplan-Meier curves for the training cohort (n = 256; log-rank P = 3.66 × 10^-6); I: Time-dependent receiver operating characteristic curves in the training cohort, with 1-, 3-, and 5-year areas under the curve of 0.792, 0.730, and 0.721, respectively; J: Kaplan-Meier curves for the internal validation cohort (n = 110; log-rank P = 0.0245); K: Time-dependent receiver operating characteristic curves in the internal validation cohort, with 1-, 3-, and 5-year areas under the curve of 0.767, 0.740, and 0.575, respectively. TCGA-LIHC: The Cancer Genome Atlas liver hepatocellular carcinoma; AUC: Area under the curve; ROC: Receiver operating characteristic.

Using the median risk score in the full TCGA-LIHC cohort, patients were divided into high- and low-risk groups, with 183 patients in each group. Kaplan-Meier analysis showed that patients in the high-risk group had significantly poorer overall survival than those in the low-risk group. Time-dependent ROC analysis in the full TCGA-LIHC cohort showed apparent discriminatory performance, with 1-, 3-, and 5-year areas under the curve (AUCs) of 0.781, 0.739, and 0.705, respectively (Figure 4C-G). Because this full-cohort analysis was based on the model-development dataset, these values were interpreted as apparent performance and were not used as the primary basis for cross-study performance comparisons.

To evaluate internal reproducibility, the TCGA-LIHC cohort was further randomly divided into a training cohort (n = 256) and an internal validation cohort (n = 110). In the training cohort, high-risk patients had significantly poorer overall survival than low-risk patients (log-rank P = 3.66 × 10^-6). The continuous risk score was significantly associated with overall survival [HR = 2.72, 95% confidence interval (CI): 1.89-3.92, P = 8.27 × 10^-8], and the 1-, 3-, and 5-year AUCs were 0.792, 0.730, and 0.721, respectively. In the internal validation cohort, the survival difference between the two risk groups remained significant (log-rank P = 0.0245). The continuous risk score was also significantly associated with overall survival (HR = 2.41, 95%CI: 1.32-4.40, P = 0.00422), with 1-, 3-, and 5-year AUCs of 0.767, 0.740, and 0.575, respectively. Harrell’s C-index was 0.713 in the training cohort and 0.654 in the internal validation cohort. These findings support the internal reproducibility of the prognostic association and indicate that discrimination was maintained at 1 and 3 years. However, the decrease in the 5-year AUC from 0.721 in the training cohort to 0.575 in the internal validation cohort suggests potential optimism in the training dataset and supports cautious interpretation of long-term predictive performance. The corresponding Kaplan-Meier and time-dependent ROC analyses are shown in Figure 4H-K.

External performance evaluation and prognostic assessment of the risk score

To further assess the external performance of the six-gene signature, analyses were performed in the independent GSE14520 and ICGC-LIRI-JP cohorts. In the GSE14520 cohort, the continuous risk score was significantly associated with overall survival in univariate Cox analysis (HR = 1.60, 95%CI: 1.01–2.54, P = 0.0474). The AUC values for 1-, 3-, and 5-year survival prediction were 0.623, 0.567, and 0.508, respectively. However, Kaplan-Meier analysis based on dichotomized risk groups did not reach statistical significance (P = 0.1577). This discrepancy may be attributed to cohort heterogeneity and the instability of dichotomized cutoff-based risk stratification, whereas the continuous risk score showed a more consistent prognostic association across datasets.

In the ICGC-LIRI-JP cohort, 93 patients with matched expression and survival data were included. The continuous risk score was also significantly associated with overall survival (HR = 2.58, 95%CI: 1.15-5.80, P = 0.0219). The 1-, 3-, and 5-year AUC values were 0.682, 0.648, and 0.581, respectively. Similarly, the dichotomized high- and low-risk groups did not show a statistically significant survival difference in the ICGC cohort (P = 0.3274). Given that ICGC-LIRI-JP represented an independent external validation cohort and that corresponding ICGC-based AUCs were available in previously reported signatures, the ICGC results were used as the main basis for descriptive cross-study performance comparisons. These results suggest that the risk score showed a consistent prognostic association across external cohorts when analyzed as a continuous variable, whereas cutoff-based discrimination remained limited because Kaplan-Meier separation did not reach statistical significance in the external cohorts (Figure 5). A detailed summary of external validation performance is provided in Table 1.

Figure 5 External performance evaluation and prognostic assessment of the six-gene risk score. A: Univariate Cox regression analysis in The Cancer Genome Atlas liver hepatocellular carcinoma; B: Multivariable Cox regression analysis adjusted for age, sex, pathological stage, and tumor grade; C: Cox analysis of the continuous risk score across The Cancer Genome Atlas liver hepatocellular carcinoma, GSE14520, and ICGC-LIRI-JP; D: Summary of time-dependent areas under the curve across datasets; E: Kaplan-Meier curve in GSE14520; F: Kaplan-Meier curve in ICGC-LIRI-JP. TCGA-LIHC: The Cancer Genome Atlas liver hepatocellular carcinoma; AUC: Area under the curve; ROC: Receiver operating characteristic; ICGC-LIRI-JP: International Cancer Genome Consortium Liver Cancer-RIKEN Japan.

Table 1 External validation performance summary of the six-gene risk score.

Dataset	n	HR	95%CI	Cox P value	KM log-rank P	1-year AUC	3-year AUC	5-year AUC
GSE14520	221	1.60	1.01-2.54	0.047	0.158	0.623	0.567	0.508
ICGC-LIRI-JP	93	2.58	1.15-5.80	0.022	0.327	0.682	0.648	0.581

Hazard ratios were derived from univariate Cox regression using the continuous risk score. The dichotomized Kaplan-Meier analysis was used only to visualize risk-group separation. AUC: Area under the curve; CI: Confidence interval; HR: Hazard ratio; KM: Kaplan-Meier; ICGC-LIRI-JP: International Cancer Genome Consortium Liver Cancer-RIKEN Japan.

In TCGA-LIHC, univariate Cox regression analysis showed that the risk score, risk group, and pathological stage were significantly associated with overall survival. In multivariable Cox regression analysis adjusted for age, sex, pathological stage, and tumor grade, the risk score remained independently associated with overall survival (HR = 2.96, 95%CI: 1.99-4.41, P = 8.08 × 10^-8). Pathological stage also remained independently associated with overall survival (HR = 2.19, 95%CI: 1.50-3.21, P = 4.90 × 10^-5) (Figure 5A and B).

A leave-one-gene-out sensitivity analysis was then performed by excluding NT5DC2. In TCGA-LIHC, the resulting five-gene score remained associated with overall survival (Kaplan-Meier log-rank P = 1.11 × 10^-7), and the 1-, 3-, and 5-year AUCs were 0.764, 0.732, and 0.700, respectively, compared with 0.781, 0.739, and 0.705 for the full six-gene score. In GSE14520, however, removal of NT5DC2 reduced the 1-, 3-, and 5-year AUCs to 0.553, 0.507, and 0.486. The continuous five-gene score was no longer significantly associated with overall survival (HR = 1.23, 95%CI: 0.67-2.26, P = 0.501), and dichotomized Kaplan-Meier separation remained nonsignificant (P = 0.144). Thus, exclusion of NT5DC2 produced only a modest reduction in apparent performance in TCGA-LIHC but a more pronounced loss of external performance in GSE14520.

Predictive performance was also placed in the context of the previously reported ferroptosis-only signature[12] and the mitochondria-only signature[13]. In accordance with the concern that AUCs derived from training or full model-development cohorts may represent optimistic apparent performance, TCGA full-cohort and training-cohort AUCs were not used as the primary basis for cross-study performance comparisons. Instead, the descriptive comparison focused mainly on validation results, particularly external ICGC validation results when corresponding data were available. In ICGC-LIRI-JP, the present six-gene model showed 1-, 3-, and 5-year AUCs of 0.682, 0.648, and 0.581, respectively. These results indicate moderate external performance, and the available external ICGC comparison did not demonstrate clear superiority over previously reported ferroptosis-only or mitochondria-only signatures. Therefore, the main contribution of the present signature lies in its compact integration of ferroptosis and mitochondrial biology and its biological interpretability, rather than in superior predictive discrimination. Because the studies used different cohort compositions, preprocessing procedures, modeling strategies, endpoint definitions, and cutoff methods, the comparison should be interpreted as descriptive context rather than as a formal statistical head-to-head comparison (Table 2).

Table 2 Descriptive comparison of external International Cancer Genome Consortium validation performance with previously reported hepatocellular carcinoma prognostic signatures.

Signature	Biological basis	Genes, n	External validation cohort	1-year AUC	3-year AUC	5-year AUC	Interpretation
Present ferroptosis-mitochondria model	Ferroptosis-mitochondria-related genes	6	ICGC-LIRI-JP	0.682	0.648	0.581	Moderate external performance
Ferroptosis-only model[12]	Ferroptosis-related genes	10	ICGC	0.680	0.718	NA	Reported external performance
Mitochondria-only model[13]	Mitochondria-related genes	6	ICGC	0.726	0.727	NA	Reported external performance

Values for References 12 and 13 were reproduced from the original reports. Training-cohort or full-cohort areas under the curve were not used as the primary basis for comparison because these estimates may represent apparent performance and may be affected by overfitting. The comparison focused on external International Cancer Genome Consortium validation results when available. Because cohort composition, expression preprocessing, endpoint definitions, modeling strategies, and cutoff methods differed across studies, this comparison should be interpreted as descriptive rather than as a formal head-to-head comparison. AUC: Area under the curve; ICGC: International Cancer Genome Consortium; NA: Not available; ICGC-LIRI-JP: International Cancer Genome Consortium Liver Cancer-RIKEN Japan.

Risk score was associated with ssGSEA-derived immune cell score patterns

To investigate associations between the prognostic signature and transcriptome-derived immune features, ssGSEA immune cell scores were compared between the high- and low-risk groups in TCGA-LIHC. After Benjamini-Hochberg correction across 25 comparisons, only plasmacytoid dendritic cell and M2 macrophage scores differed significantly between the groups; both were lower in the high-risk group (FDR = 0.0010 and 0.0199, respectively) (Figure 6A).

Figure 6 Single-sample gene set enrichment analysis-derived immune cell score patterns associated with the risk score. A: Boxplots with individual patient observations for the two immune cell scores that differed significantly between the high- and low-risk groups after Benjamini-Hochberg correction; B: Spearman correlations between the continuous risk score and immune cell scores that remained significant after false discovery rate correction; C: Sample-level heatmap of the five immune cell scores that remained significant after false discovery rate correction, ordered by risk group and continuous risk score. FDR: False discovery rate; ssGSEA: Single-sample gene set enrichment analysis.

Spearman analysis identified five immune cell scores that remained significantly correlated with the continuous risk score after FDR correction. The risk score was negatively correlated with the plasmacytoid dendritic cell score (rho = -0.231, FDR = 0.000217), the M2 macrophage score (rho = -0.157, FDR = 0.0240), the natural killer cell score (rho = -0.156, FDR = 0.0240), and the plasma cell score (rho = -0.137, FDR = 0.0440), and positively correlated with the eosinophil score (rho = 0.150, FDR = 0.0248). The modest correlation magnitudes suggest immune score remodeling rather than providing direct evidence of large changes in immune cell abundance (Figure 6B and C).

Pathway enrichment analysis revealed proliferative and metabolic features of the high-risk group

To explore the biological mechanisms underlying the prognostic risk signature, Hallmark GSVA and GSEA analyses were performed between the high- and low-risk groups. GSVA showed that the high-risk group had significantly greater activity in pathways related to DNA repair, E2F targets, the G2/M checkpoint, MYC targets V1, MYC targets V2, the unfolded protein response, the mitotic spindle, phosphoinositide 3-kinase (PI3K)-protein kinase B (AKT)-mechanistic target of rapamycin (mTOR) signaling, mTORC1 signaling, and protein secretion. In contrast, the low-risk group showed greater enrichment in bile acid metabolism, xenobiotic metabolism, fatty acid metabolism, coagulation, KRAS signaling down, complement, interferon-gamma response, inflammatory response, TNF-α signaling via nuclear factor-κB, and IL-6-JAK-STAT3 signaling (Figure 7A and B).

Figure 7 Hallmark pathway enrichment analysis of the ferroptosis-mitochondria-related risk signature. A: Gene set variation analysis heatmap of differentially enriched Hallmark pathways; B: Boxplots of representative Hallmark pathways; C: Bar plot of normalized enrichment scores from gene set enrichment analysis; D-F: Gene set enrichment analysis plots for MYC targets V1, E2F targets, and MYC targets V2. GSEA: Gene set enrichment analysis; GSVA: Gene set variation analysis.

Consistent with the GSVA findings, Hallmark GSEA similarly showed significant enrichment of proliferative and cell cycle-related pathways in the high-risk group. The most significantly enriched pathways included MYC targets V1 (NES = 2.204, adjusted P = 7.67 × 10^-10), E2F targets (NES = 2.161, adjusted P = 7.67 × 10^-10), MYC targets V2 (NES = 2.002, adjusted P = 7.67 × 10^-10), G2/M checkpoint (NES = 1.995, adjusted P = 7.67 × 10^-10), DNA repair (NES = 1.711, adjusted P = 7.67 × 10^-10), and mitotic spindle (NES = 1.634, adjusted P = 7.67 × 10^-10). Oxidative phosphorylation, PI3K-AKT-mTOR signaling, mTORC1 signaling, unfolded protein response, protein secretion, and glycolysis were also significantly enriched in the high-risk group (Figure 7C-F).

Together, the GSVA and GSEA findings suggest that the high-risk phenotype is characterized by increased proliferative activity, DNA repair, cell cycle progression, mitochondrial energy metabolism, and stress-response programs, whereas the low-risk phenotype is characterized by greater activity of liver-associated metabolic programs.

Nomogram construction, calibration, and decision curve analysis

To further explore the potential clinical predictive value of the six-gene risk score, a prognostic nomogram was constructed by integrating the risk score with age, sex, pathological stage, and tumor grade. After patients with missing clinical information were excluded, 340 patients with 114 deaths were included in the final full-cohort analysis. In the multivariable Cox model, the risk score remained strongly associated with overall survival (HR = 2.96, 95%CI: 1.99-4.41, P = 8.08 × 10^-8), and pathological stage also remained an independent prognostic factor (HR = 2.19, 95%CI: 1.50-3.21, P = 4.90 × 10^-5). Age, sex, and tumor grade were not statistically significant after adjustment but were retained a priori for clinical adjustment and comparability.

A single multivariable Cox model was used to generate separate 1-, 3-, and 5-year overall survival probability scales in the nomogram (Figure 8A). Bootstrap calibration with 200 resamples showed the relationship between model-predicted and observed survival probabilities at the three time points (Figure 8B-D). Decision curve analysis compared the potential clinical net benefit of the risk score-only, clinical, stage-only, and combined risk-plus-clinical models (Figure 8E-G).

Figure 8 Nomogram construction and validation for the six-gene risk score in The Cancer Genome Atlas liver hepatocellular carcinoma. A: Single multivariable Cox nomogram integrating risk score, age, sex, pathological stage, and tumor grade, with separate scales for 1-, 3-, and 5-year overall survival probabilities; B-D: Bootstrap calibration curves for 1-, 3-, and 5-year overall survival prediction; E-G: Decision curve analyses at 1, 3, and 5 years; H: Harrell’s C-index for the risk score-only and combined models in the training and internal validation cohorts.

Holdout internal validation was performed using the predefined TCGA split after complete-case filtering (training cohort, n = 238; internal validation cohort, n = 102). In the training cohort, Harrell’s C-index was 0.708 for the risk score-only model and 0.742 for the combined nomogram model. In the internal validation cohort, the corresponding C-index values were 0.662 and 0.669, respectively (Figure 8H). The attenuation in the holdout cohort indicates moderate rather than definitive generalizability, whereas the small increase from 0.662 to 0.669 suggests only limited incremental discrimination from the clinical variables. After this validation, the nomogram shown in Figure 8A was refitted using all 340 complete-case TCGA patients for presentation.

Exploratory protein-level assessment of the six-gene signature using public CPTAC data

To provide complementary protein-level evidence for the six-gene signature, public CPTAC proteomic data accessed through UALCAN were used to compare protein abundance between HCC and normal liver tissues. Among the six model genes, NT5DC2, TIMM9, and LIG3 showed significantly higher protein abundance in HCC tissues than in normal liver tissues. In contrast, LONP1 and PRDX6 showed significantly lower protein abundance in HCC tissues, whereas no significant difference was observed for CISD1.

Specifically, NT5DC2 (P = 3.92 × 10^-31), TIMM9 (P = 1.21 × 10^-31), and LIG3 (P = 1.64 × 10^-74) showed significantly higher protein abundance in HCC tissues, whereas LONP1 (P = 1.41 × 10^-4) and PRDX6 (P = 3.73 × 10^-45) showed significantly lower protein abundance. CISD1 protein abundance did not differ significantly between HCC tissues and normal liver tissues (P = 0.921). These findings showed heterogeneous and partially discordant protein expression patterns across the six genes in HCC (Supplementary Figure 1; Supplementary Table 1).

DISCUSSION

In this study, we systematically investigated ferroptosis-mitochondria-related genes in HCC and evaluated their relevance to molecular classification, prognostic stratification, immune characterization, and pathway alterations. By intersecting ferroptosis-related genes from FerrDb with mitochondrial genes from MitoCarta 3.0, we obtained 47 ferroptosis-mitochondria-related genes. Based on these genes, four molecular subtypes with distinct survival outcomes were identified in TCGA-LIHC, and a six-gene prognostic signature comprising NT5DC2, TIMM9, CISD1, LONP1, PRDX6, and LIG3 was constructed. The risk score was independently associated with overall survival in TCGA-LIHC and was also significantly associated with overall survival as a continuous variable in both GSE14520 and ICGC-LIRI-JP. These findings suggest that the integrated ferroptosis-mitochondria-related gene set may represent a biologically interpretable aspect of HCC heterogeneity that links metabolic stress, the tumor immune microenvironment, pathway dysregulation, and prognostic risk.

A key feature of the present study is that ferroptosis and mitochondrial dysfunction were analyzed jointly rather than in isolation. Ferroptosis is closely related to iron metabolism, lipid peroxidation, redox imbalance, and mitochondrial activity, whereas mitochondrial dysfunction is involved in energy metabolism, reactive oxygen species generation, apoptotic signaling, and metabolic reprogramming. The intersection of ferroptosis-related and mitochondrial genes therefore yielded a focused gene set for investigating HCC heterogeneity in the context of metabolic stress and regulated cell death. The four molecular subtypes identified based on these 47 genes indicate that ferroptosis-mitochondria-related transcriptional patterns may reflect clinically relevant molecular diversity in HCC.

Molecular subtype analysis further indicated that ferroptosis-mitochondria-related genes were associated with both prognosis and immune heterogeneity. The four subtypes showed significantly different overall survival, and broad differences in ssGSEA-derived immune cell scores were observed among the four subtypes. In particular, C3 showed the lowest scores for many immune cell populations and therefore exhibited an immune-low score pattern. This pattern may be related to reduced transcriptome-derived immune activity rather than changes in a single immune lineage. It may also be consistent with reported ferroptosis–immune crosstalk, including CD8+ T-cell-mediated induction of tumor ferroptosis and ferroptosis-associated immune evasion[6,7]. The molecular subtype results therefore provide a basis for future comparison with established immune-excluded or immune-cold HCC classifications and for investigating their relationship with ferroptosis-mitochondria-related stress programs. This is consistent with previous evidence that HCC can be stratified into molecularly and immunologically distinct subgroups[35,36].

The six-gene risk signature showed prognostic relevance and moderate predictive performance in TCGA-LIHC. High-risk patients had significantly poorer overall survival, and the risk score remained independently associated with survival after adjustment for age, sex, stage, and grade. In the external GSE14520 and ICGC-LIRI-JP cohorts, the continuous risk score was significantly associated with overall survival, supporting a consistent trend across datasets. However, dichotomized high- and low-risk stratification based on the selected cutoff did not reach statistical significance in these external cohorts. This discrepancy may be partly attributable to cross-platform expression differences, cohort-specific clinical composition, limited sample size, and the instability of cohort-specific median cutoffs used only for Kaplan-Meier visualization. Therefore, the signature should be interpreted as a prognosis-associated score rather than as an already optimized clinical classifier. The corrected concordance analysis yielded C-index values of 0.713 and 0.654 in the training and internal validation cohorts, respectively, indicating moderate discrimination with attenuation in the holdout subset. This pattern was consistent with the lower 5-year AUC observed in the internal validation cohort and supports cautious interpretation of long-term predictive performance. Importantly, AUCs derived from the training or full model-development cohort may overestimate generalizable performance because of possible overfitting. This concern was reflected by the decrease in the 5-year AUC from 0.721 in the training cohort to 0.575 in the internal validation cohort. Therefore, the training and full-cohort TCGA results should be regarded as apparent performance, whereas the internal and external validation results provide a more appropriate basis for evaluating model generalizability.

The genes included in the signature also have biologically plausible links to tumor progression and mitochondrial stress. TIMM9 and LONP1 are involved in mitochondrial protein transport[37] and mitochondrial proteostasis[38], respectively. CISD1 and PRDX6 have been implicated in mitochondrial redox regulation, iron-sulfur homeostasis, antioxidant defense, and lipid peroxidation-related stress, all of which may influence ferroptosis sensitivity[39,40]. LIG3 is linked to DNA repair capacity and mtDNA integrity[41], which is consistent with the enrichment of DNA repair and proliferative pathways in high-risk tumors. NT5DC2 has been implicated in nucleotide metabolism and HCC progression and was previously reported to promote tumor cell proliferation by stabilizing EGFR[42]. Although these mechanisms were not experimentally tested in the present study, the composition of the signature suggests that multiple layers of mitochondrial homeostasis, redox balance, and stress adaptation may contribute to HCC progression.

The NT5DC2 ablation analysis provided additional context for its contribution to the signature. Removing NT5DC2 caused only a modest decrease in TCGA AUCs but a clearer reduction in GSE14520 discrimination, where the continuous five-gene score lost statistical significance and the 5-year AUC fell below 0.50. This pattern supports a meaningful contribution of NT5DC2 to the transportability of the score. However, because the ablation analysis retained the other coefficients without refitting and was evaluated in a cross-platform external dataset, it does not establish NT5DC2 as the sole driver of model performance or prove a causal role.

Direct descriptive comparison with the previously reported ferroptosis-only[12] and mitochondria-only[13] signatures did not demonstrate superior predictive performance for the present model. In response to the concern that training-cohort or full model-development cohort AUCs may overestimate generalizable performance, we interpreted TCGA full-cohort and training-cohort results as apparent performance and did not use them as the main basis for cross-study comparison. Based on external ICGC validation, the present model showed moderate performance, with 1-, 3-, and 5-year AUCs of 0.682, 0.648, and 0.581, respectively. The available 1- and 3-year ICGC AUCs were broadly comparable to, but not consistently higher than, those reported for ferroptosis-only and mitochondria-only signatures. Therefore, the added value of the present signature lies primarily in its compact integration of ferroptosis and mitochondrial biology, its biological interpretability, and its multi-cohort evaluation rather than in a claim of superior discrimination. Because these AUCs were obtained using different preprocessing and modeling pipelines, they should not be interpreted as a formal statistical head-to-head comparison. Other single-domain HCC signatures, including pyroptosis- and inflammatory-response-related models, further illustrate the diversity of prognostic modeling strategies[43,44].

Immune cell score analysis further supported the biological relevance of the risk score. The risk score was negatively correlated with scores for plasmacytoid dendritic cells, M2 macrophages, natural killer cells, and plasma cells, and positively correlated with the eosinophil score. Although the absolute correlation coefficients were modest, these findings suggest that the high-risk phenotype may be accompanied by immune microenvironment remodeling. The strongest negative correlation was observed for the plasmacytoid dendritic cell score, suggesting a potential relationship between ferroptosis-mitochondria-related risk status and antigen-presenting or interferon-related immune components, possibly through ferroptosis-driven release of signaling molecules that remodel the tumor microenvironment[6,7]. These computational findings should be interpreted cautiously, because ssGSEA estimates relative immune cell activity indirectly from bulk transcriptomic profiles and cannot replace direct immune profiling.

Functional enrichment analyses provided a coherent biological interpretation of the high-risk phenotype. Hallmark GSVA and GSEA consistently showed that high-risk tumors were enriched in MYC targets, E2F targets, G2/M checkpoint, DNA repair, the mitotic spindle, the unfolded protein response, PI3K-AKT-mTOR signaling, mTORC1 signaling, oxidative phosphorylation, and glycolysis. These pathways are involved in cell cycle progression, transcriptional activation, proliferative capacity, genomic maintenance, mitochondrial energy metabolism, and cellular stress responses. The enrichment of both oxidative phosphorylation and glycolysis suggests metabolic plasticity rather than a simple shift toward one energy pathway. Together with MYC/E2F activation and mTOR signaling, this pattern may reflect a tumor state characterized by increased proliferative demand, adaptation to mitochondrial stress, and altered susceptibility to oxidative injury and ferroptotic stress. In contrast, low-risk tumors showed greater enrichment in bile acid metabolism, xenobiotic metabolism, fatty acid metabolism, coagulation, complement, and inflammatory response pathways, suggesting that the low-risk phenotype may retain stronger liver-associated metabolic programs.

To explore the potential clinical predictive value, we constructed a multivariable Cox nomogram with separate 1-, 3-, and 5-year survival probability scales. Age, sex, and tumor grade were retained a priori to preserve clinical adjustment and comparability, although they were not independently significant in the full model. Bootstrap calibration assessed optimism-corrected agreement, and holdout internal validation showed that the combined model C-index declined from 0.742 in the training cohort to 0.669 in the internal validation cohort. The clinical variables produced only a small incremental improvement over the risk score alone in the holdout data (0.669 vs 0.662). These results support the nomogram as a preliminary prediction framework but also indicate moderate discrimination and the need for independent external calibration and prospective validation before clinical use.

Public CPTAC-based exploration provided complementary protein-level observations, although the patterns were not uniformly concordant with the transcriptomic findings. NT5DC2, TIMM9, and LIG3 were significantly increased in HCC tissues, whereas LONP1 and PRDX6 were decreased; CISD1 showed no significant difference. These findings suggest that several model genes may be dysregulated at the protein level in HCC. Nevertheless, the protein expression patterns did not fully align with the prognostic directions of all genes, and the public CPTAC-based analysis should be regarded as exploratory rather than as experimental validation. The observed transcript-protein discordance may reflect post-transcriptional regulation, differences in translational efficiency, protein turnover or degradation, and intratumoral cellular heterogeneity. Moreover, tumor-versus-normal differential protein expression and the association of interpatient transcriptomic variation with survival represent distinct biological questions; therefore, a nonsignificant or directionally discordant tumor-normal protein difference does not necessarily invalidate the prognostic association identified at the transcriptomic level. Further in-house protein assays and functional experiments are needed to determine whether these genes directly regulate ferroptosis, mitochondrial dysfunction, and HCC progression.

Several limitations should be acknowledged. First, this study was based mainly on retrospective public transcriptomic datasets. Gene-wise z-standardization within each cohort reduced cross-platform scale differences, but because the cohorts were analyzed separately and were not pooled for ComBat correction, residual microarray/RNA-sequencing batch effects may still have contributed to the weaker external performance. Etiological heterogeneity may also be relevant because TCGA contains a mixed population, GSE14520 is predominantly derived from Chinese patients with hepatitis B virus-related HCC, and ICGC-LIRI-JP represents a Japanese cohort with different proportions of hepatitis B virus-, hepatitis C virus-, and nonviral HCC. Second, although the continuous risk score remained significantly associated with overall survival in both external cohorts, dichotomized high- and low-risk groups did not show significant survival differences, and the GSE14520 5-year AUC of 0.508 was close to chance-level discrimination. The 7:3 holdout analysis reduced, but did not eliminate, concerns about overfitting. In particular, the 5-year AUC decreased from 0.721 in the training cohort to 0.575 in the internal validation cohort, indicating that long-term predictive performance was less stable than short- and intermediate-term performance. The nomogram C-index also declined from training to validation. Therefore, the model should not be considered a definitive long-term prognostic classifier without further independent and prospective validation. Third, the NT5DC2 ablation analysis was a coefficient-preserving sensitivity analysis rather than a fully refitted competing model, and its findings should be interpreted as evidence of contribution rather than causality. Although pathway enrichment was performed using genome-wide transcriptomic profiles rather than only the six model genes, risk-group assignment and downstream transcriptional states were not entirely independent. Accordingly, the GSVA and GSEA findings should be interpreted as associative and hypothesis-generating rather than as causal evidence that the risk score directly drives the observed pathway alterations. Fourth, public CPTAC analysis provided only exploratory protein-level observations, and no in-house experimental validation was performed. Treatment and recurrence information was incomplete across the public datasets, preventing evaluation of treatment response prediction. Finally, immune cell and pathway scores were inferred computationally from bulk transcriptomic data and require validation using direct immune profiling, single-cell approaches, and prospective clinical cohorts.

CONCLUSION

This study identified ferroptosis-mitochondria-related molecular subtypes and developed a six-gene prognostic signature in HCC. The risk score was associated with overall survival, immune cell score patterns, proliferative and metabolic pathway alterations, and heterogeneous protein-level expression patterns in public CPTAC data. The validation results suggest moderate rather than definitive predictive performance. These findings indicate that ferroptosis-mitochondria-related genes may provide useful clues for prognostic assessment and biological interpretation. Further experimental and prospective clinical validation is required before clinical application.