Gene Signature Tied to Folate Metabolism Predicts Liver Cancer Survival Unevenly

Researchers in China have built a 20-gene prognostic model for hepatocellular carcinoma (HCC), the most common form of primary liver cancer, drawing on genes linked to folate metabolism. The signature, developed from transcriptomic data in The Cancer Genome Atlas liver cohort (TCGA-LIHC), showed strong apparent discrimination in the training set but a more modest performance after statistical correction, and its results varied sharply across two independent validation cohorts. The work, published in Frontiers in Genetics, illustrates both the promise and the pitfalls of computational prognostic signatures for a disease that remains notoriously difficult to stage at the molecular level.

Hepatocellular carcinoma accounts for more than 90 percent of primary liver cancers and carries a five-year survival rate below 20 percent, according to the epidemiological figures cited by the study team. Current treatment options range from surgical resection and local ablation to transarterial chemoembolization, systemic targeted therapy, and immunotherapy. Yet the standard anatomical staging systems, such as TNM and BCLC, capture little of the disease’s molecular heterogeneity, and patients at similar clinical stages can follow markedly different courses. This gap has driven a wave of efforts to build transcriptome-based risk models that might complement clinical indicators, and the new study adds folate metabolism to the list of biological angles explored so far, which has included immune-related genes and broadly defined metabolic gene sets.

Folate metabolism sits at the heart of the one-carbon metabolic network. By supplying one-carbon units, it feeds purine and thymidylate synthesis, generates methyl donors, and helps maintain cellular redox balance, all processes that support proliferation. Disruptions in folate-cycle enzymes and one-carbon metabolism have been implicated in tumor initiation and progression in the liver, and epidemiological studies have linked lower folate levels to liver injury and increased HCC risk. To probe this connection systematically, the researchers mined the GeneCards database for protein-coding genes associated with folate metabolism, retrieving 2,043 candidates with relevance scores above 1.

Using the TCGA-LIHC cohort, which comprised 371 tumor samples and 50 normal controls, the team performed differential expression analysis with the limma empirical Bayes method and Benjamini-Hochberg false discovery rate control. The thresholds were an adjusted p-value below 0.05 and an absolute log2 fold change above 0.5. The screen yielded 974 differentially expressed folate metabolism-related genes, 753 upregulated and 221 downregulated in tumors. Enrichment analysis showed these genes clustered around DNA replication, DNA recombination, double-strand break repair, telomere maintenance, and responses to oxidative stress and hypoxia, alongside pathways including the cell cycle, HIF-1 and p53 signaling, homologous recombination, and, notably, one-carbon metabolism by folate and folate transport itself.

For model construction, the researchers matched gene expression with complete survival data for 318 HCC patients, among whom 108 death events occurred. Univariate Cox regression with false discovery correction identified 231 genes significantly associated with overall survival. Elastic Net Cox regression, with the mixing parameter alpha set at 0.9 and the penalty lambda minimized at 0.056 through 10-fold cross-validation, then whittled the field down to 20 genes. The resulting risk score weights EPO, MYCN, PPAT, PRDX1, MAPT, TALDO1, GPX4, SLC2A1, SPP1, EZH2, PSMA7, DTYMK, LGALS3, HRAS, BMI1, BMP2, and PGF positively, while PON1, LCAT, and SERPINF2 carry negative coefficients. A LASSO Cox model built on the same candidate set selected the identical 20 genes and achieved a nearly identical C-index, suggesting the signature is robust to the choice of penalty method.

In the training cohort, the model’s Harrell concordance index reached 0.767, with a 95 percent confidence interval of 0.718 to 0.815. Time-dependent ROC analysis gave area-under-the-curve values of 0.862, 0.768, and 0.737 for one-, three-, and five-year overall survival, respectively, and at five years the risk score outperformed clinical stage, age, and sex. Multivariate Cox regression adjusting for age, sex, and stage confirmed the risk score as an independent prognostic factor, with a hazard ratio of 4.14. But the authors were candid about overfitting risk: after 999 rounds of bootstrap resampling in which the entire modeling procedure was repeated, the optimism-corrected C-index fell to 0.676, and with only 108 events spread across 20 predictors, the events-per-variable ratio of 5.4 left little margin. The 231 candidate genes preceding penalization corresponded to just 0.47 events per variable.

To interpret the model, the team applied SHAP analysis, computing each gene’s contribution as its regression coefficient times the deviation of a patient’s expression from the cohort mean. EPO, MYCN, PPAT, PON1, and PRDX1 showed the highest global importance, with mean absolute SHAP values of 0.173, 0.142, 0.105, 0.103, and 0.083. Elevated expression of EPO, MYCN, and PRDX1 generally pushed risk scores upward, while higher PON1 expression pulled them down, consistent with its negative coefficient. For the single patient with the highest risk score in the cohort, MYCN contributed the largest individual SHAP value at +0.696, followed by EPO, PON1, and PPAT.

External validation told a more complicated story. In the LIRI-JP cohort of 68 Japanese patients sequenced by RNA-seq, the model achieved a C-index of 0.735 and separated high- from low-risk groups with a log-rank p-value of 0.0075. In GSE10141, however, a cohort of 62 evaluable patients profiled on a DASL 6K platform using formalin-fixed paraffin-embedded samples, the C-index dropped to 0.495, essentially chance, and the survival difference did not reach significance. The authors note that GSE10141 recorded no deaths during the first year of follow-up and that platform and sample-processing differences may underlie the discrepancy, though the exact causes remain unexplained. Risk stratification thresholds and regression coefficients were held fixed from the training cohort, with no re-optimization in the validation sets.

Beyond survival prediction, the study mapped immune and drug-response features onto the risk groups using computational tools. CIBERSORT deconvolution estimated that eight of 22 immune cell types differed in relative proportion between risk groups, with higher estimated fractions of M0 macrophages, regulatory T cells, and follicular helper T cells in the high-risk group. Immunophenoscore values were significantly higher in the low-risk group, while TIDE analysis showed lower T-cell dysfunction but higher T-cell exclusion scores in high-risk patients. Crucially, the proportion of patients predicted by TIDE to respond to immune checkpoint inhibitors did not differ significantly between groups, at 84.9 percent versus 82.4 percent, so the authors caution that these immune-score differences cannot be read as differences in actual immunotherapy benefit. OncoPredict analysis based on GDSC2 data flagged 123 of 198 drugs with differing predicted sensitivities, though these remain computational estimates requiring pharmacological confirmation.

A parallel line of work focused on GSS, the gene encoding glutathione synthetase, which screened as differentially expressed and prognostic but was not selected for the final model; adding it back raised the C-index only from 0.767 to 0.768 without statistical significance. Nevertheless, GSS expression was elevated in HCC tissues and high expression correlated with shorter survival in both TCGA-LIHC and LIRI-JP. Single-cell analysis of the GSE149614 dataset identified malignant hepatocytes by inferred copy-number variation and found GSS expression markedly higher in malignant cells, with the fraction of GSS-positive cells rising from 27.1 to 44.1 percent. In vitro, knocking down GSS in Huh-7 cells impaired proliferation, migration, and invasion, while overexpression in SNU-182 cells produced the opposite effects. The authors emphasize that they did not measure glutathione levels, reactive oxygen species, or metabolic flux, so the mechanism remains open. Taken together, the findings support GSS as a candidate driver of malignant behavior and underscore that the 20-gene signature, however promising in RNA-seq cohorts, still needs multicenter validation before any clinical use.

Subject of Research: A transcriptomic prognostic model based on folate metabolism-related genes for hepatocellular carcinoma survival prediction and functional analysis of GSS

Article Title: A folate metabolism-related prognostic signature for predicting survival in patients with hepatocellular carcinoma

Article References: Lin, L., Chen, X., Wu, J., Li, C., Yu, M., & Cai, B. (2026). A folate metabolism-related prognostic signature for predicting survival in patients with hepatocellular carcinoma. Frontiers in Genetics, 17, Article 1822019. https://doi.org/10.3389/fgene.2026.1822019

Image Credits: AI Generated

DOI: 10.3389/fgene.2026.1822019

Keywords: hepatocellular carcinoma, folate metabolism, prognostic signature, TCGA-LIHC, Elastic Net Cox regression, SHAP analysis, GSS, glutathione synthetase, single-cell RNA sequencing, tumor immune microenvironment, drug sensitivity prediction, biomarker validation