Characterization of an Activated Metabolic Transcriptional Program in Hepatoblastoma Tumor Cells Using scRNA-seq

Abstract

Hepatoblastoma is the most common primary liver malignancy in children, with metabolic reprogramming playing a critical role in its progression due to the liver's intrinsic metabolic functions. Enhanced glycolysis, glutaminolysis, and fatty acid synthesis have been implicated in hepatoblastoma cell proliferation and survival. In this study, we screened for altered overexpression of metabolic enzymes in hepatoblastoma tumors at tissue and single-cell levels, establishing and validating a hepatoblastoma tumor expression metabolic score using machine learning. Starting from the Mammalian Metabolic Enzyme Database, bulk RNA sequencing data from GSE104766 and GSE131329 datasets were analyzed using supervised methods to compare tumors versus adjacent liver tissue. Differential expression analysis identified 287 significantly regulated enzymes, 59 of which were overexpressed in tumors. Functional enrichment in the KEGG metabolic database highlighted a network enriched in amino acid metabolism, as well as carbohydrate, steroid, one-carbon, purine, and glycosaminoglycan metabolism pathways. A metabolic score based on these enzymes was validated in an independent cohort (GSE131329) and applied to single-cell transcriptomic data (GSE180665), predicting tumor cell status with an AUC of 0.98 (sensitivity 0.93, specificity 0.94). Elasticnet model tuning on individual marker expression revealed top tumor predictive markers, including FKBP10, ATP1A2, NT5DC2, UGT3A2, PYCR1, CKB, GPX7, DNMT3B, GSTP1, and OXCT1. These findings indicate that an activated metabolic transcriptional program, potentially influencing epigenetic functions, is observed in hepatoblastoma tumors and confirmed at the single-cell level.

Keywords: DNA methylation; cancer; epigenetics; hepatoblastoma; metabolism; one-carbon.

Figure 1

Differentially expressed enzymes in tumors compared to normal adjacent liver tissues from human hepatoblastoma samples: (A) Volcano plot of differentially expressed enzymes in the GSE104766 transcriptome dataset; (B) Principal component analysis of 287 differentially expressed enzymes; (C) Unsupervised clustering (Euclidean distances) with an expression heatmap of the 287 differentially expressed enzymes.

Figure 2

Activated metabolic transcriptional program in hepatoblastoma tumors: GSE104766: (A) Bar plot of functional enrichment (KEGG database) performed on 59 enzymes found overexpressed in hepatoblastoma tumors; (B) Emaplot of functional enrichment (KEGG database) performed on 59 enzymes found overexpressed in hepatoblastoma tumors.

Figure 3

Activated metabolic network in human hepatoblastoma tumors: GSE104766, Functional enrichment network of metabolic enzyme actors found over expressed in human hepatoblastoma tumors (enrichment on KEGG database), color scale from yellow to purple represent the negative logarithm base 10 of the enriched p-values.

Figure 4

Validation of the metabolic signature on an external cohort of hepatoblastoma transcriptomes: dataset GSE131329. (A) Unsupervised clustering and expression heatmap of the metabolic score signature. (B) Calculation of the metabolic score according to tissue groups: the p-value evaluating the difference between the two tissue groups was obtained by a two-sided Student’s t-test.

Figure 5

High level of metabolic score at the single-cell level in tumor cells from hepatoblastoma: GSE180665 (A) UMAP dimensional reduction with expression of Glypican 3 (GPC3). (B) Violin plot of the metabolic score across sample groups (p-value, two-sided t-test). (C) Feature plot of metabolic score, split by sample types (computed on 59 HB-activated enzymes) on UMAP dimensional reduction. (D) Violin plot of metabolic score stratified by cell types: normal hepatocytes (Hep) in pink. (E) Violin plot of metabolic score comparing normal hepatocytes in blue (Hep) versus tumor cells (PDX + human tumors) by scRNA-seq. (F) ROC curve and area under the curve for metabolic score in predicting hepatoblastoma tumor cells as compared to normal hepatocytes.

Figure 6

Higher levels of metabolic score quantified in G1 quiescent cells derived from human hepatoblastoma tumoroids: scRNA-seq dataset GSE233923. (A) t-SNE dimensional reduction of scRNA-seq transcriptomes from five human hepatoblastoma tumoroids post canonical correlation integration. (B) t-SNE dimensional reduction with cell cycle phase prediction (G1: Gap1, S: DNA synthesis, G2M: Gap2/Mitosis). (C) Metabolic score quantification on t-SNE dimensional reduction with stratification by cell cycle phases (p-values in tables were obtained by a two-sided Student’s t-test in the table). (D) Dot plot stratified by cell cycle phases for quantification of metabolic score, KRT19, and GPC3 expression.

Figure 7

Importance of overexpressed metabolic markers to predict tumor cell status in scRNA-seq of hepatoblastoma samples: GSE180665. (A) Alpha and lambda parameter tuning from the ElasticNet model with a 70/30 training/validation split of the GSE180665 dataset: tumor cells versus normal hepatocyte status prediction. (B) Best lambda estimation for the alpha parameter fixed at 0.1 in the ElasticNet model. (C) Positive ElasticNet coefficients for 41 metabolic markers expressed in HB tumor cells at the single-cell level, with individual area under the curve (AUC) values for the ten best predictive enzymes.