Comprehensive characterization of β-alanine metabolism-related genes in HCC identified a novel prognostic signature related to clinical outcomes

Data collection and processing

To conduct this study, we obtained gene expression data and clinical annotations from two well-known databases: GEO and TCGA. To ensure the reliability of our analysis, we excluded patients without survival information or with survival times less than 30 days. Additionally, we used the ComBat method from the “SVA” package to address any batch effects caused by nonbiological technical biases [13]. The dataset for our study consisted of 569 HCC patients, incorporating information from TCGA-LIHC, GSE76427, GSE116174, and GSE144269 datasets. This comprehensive collection of expression profiles and clinical data enabled us to conduct a thorough and in-depth investigation of HCC. From the KEGG_BETA_ALANINE_METABOLISM dataset in the MsigDb platform, we extracted 22 classic βAMRGs.

The single-cell RNA (scRNA) sequencing data for HCC was sourced from the GSE166635 dataset [14], which included two HCC single-cell samples (GSM5076749 and GSM5076750). The Read10X function was employed to read the single-cell data, with the following filtering criteria set: min.cells = 3, min.features = 200, nCount_RNA ≥ 1000, nFeature_RNA ≥ 200, nFeature_RNA ≤ 9000, percent.mt ≤ 20.

Pan-cancer analysis

To gain insights into the molecular features of βAMRGs in different human cancers, we conducted a thorough pan-cancer analysis, integrating genomic, transcriptomic, and clinical data [15, 16]. Our initial step involved gathering and consolidating raw data and clinical information from pan-cancer cohorts via the Xena website. The CNV and SNV data from the TCGA database were then processed and visualized as heatmaps, offering a comprehensive overview of variations in βAMRGs across various cancers. Furthermore, we conducted pan-cancer analyses on methylation levels and differential mRNA expressions.

To evaluate the prognostic significance of βAMRGs in diverse malignancies, we performed univariate COX regression analysis on their expression in relation to OS. By applying these extensive methodologies, our objective was to uncover the possible functions and importance of β-alanine metabolism-related genes in various cancer types, thus providing new perspectives and potential therapeutic avenues in the realm of cancer research.

scRNA sequencing analysis

After single-cell quality control, we initially employed the FeatureScatter function to visualize the quality control results. The LogNormalize method was applied for data standardization [17], the FindVariableFeatures function for downstream dimensionality reduction analysis, the ScaleData function for data normalization, and the RunHarmony function for batch correction and dimensionality reduction analysis. The JackStraw function was utilized to determine the appropriate PC values, while the singleR function facilitated automated cell annotation. The AUCell, Add, singscore, ssGSEA, and UCell algorithms were employed to assess the metabolic states of each cell, and the cumulative results of all algorithm scores were referred to as ‘total scoring’ in this study, abbreviated as Scoring.

Violin plots were used to illustrate the differences in metabolic activity between each cell type. Finally, the CopyKat algorithm was employed to predict the benign or malignant nature of each cell. The Wilcoxon test function was used to compare metabolic differences between benign and malignant cells.

Clustering analysis

Due to significant variations in gene expression profiles across the collected datasets, we devised a β-alanine metabolism score model to emphasize the distinct expression levels among samples. Initially, we employed the ssGSEA algorithm to calculate β-alanine metabolism enrichment scores for all 569 patients, thus assigning each patient a specific β-alanine metabolism score indicative of pathway activity.

Next, we conducted differential analysis, enabling us to explore expression pattern differences among the samples. These results were visually represented in a heatmap, generated using the “pheatmap” package, providing a comprehensive illustration of the clustering analysis outcomes.

By comparing mRNA expression levels of genes in tumor tissues with those in normal tissues, we successfully categorized mRNA expression statuses in tumor tissues into three distinct groups. This approach provided valuable insights into the dynamic nature of β-alanine metabolism across various cancer types, potentially revealing implications for cancer research and therapeutic strategies.

Prognostic performances and molecular characteristics of β-alanine metabolism subtypes

To delve deeper into the connections between gene expression levels within these three clusters, we skillfully created violin plots using the “ggpubr” package to illustrate the enrichment scores among them. To examine prognostic differences among these clusters, we employed the “survival” and “survminer” packages in R.

By leveraging the resources of the MsigDb platform and drawing insights from previous research, we pinpointed 42 conventional metabolic pathways, 24 immune-related pathways, and 10 cell death pathways. To gauge the activity of these pathways within each subtype of β-alanine metabolism, we employed the GSVA program to calculate the metabolic score, immunological score, and cell death score for each liver cancer sample.

To assess variations in pathway activity among the three subtypes, we carried out the Kruskal-Wallis test. Furthermore, a comprehensive investigation was conducted to explore alterations in immune cell infiltration and expression of immune checkpoint genes (ICGs), providing detailed insights into the changes within the immune environment across diverse subtypes. For this analysis, we utilized various immunological algorithms available on the TIMER2.0 platform, including TIMER, CIBERSOFT, QUANTISEQ, EPIC, and others [18, 19]. Moreover, we applied the Kruskal-Wallis test to investigate variations in immune cell infiltration and expression of ICGs across the subtypes. We exclusively showcased the outcomes that demonstrated statistical significance with a P-value lower than 0.05.

β-alanine metabolism is an auxiliary clinical indicator to predict the prognosis of HCC patients

The TCGA dataset contained 343 HCC samples, which were randomly divided into two groups: the training cohort, comprising 60% of the individuals, and the test1 cohort, including 40% of the participants. Moreover, all samples from the TCGA dataset were assigned to the test2 cohort (n=343), while the test3 cohort (n=226) exclusively comprised samples from the GEO platform.

To address collinearity and prevent overfitting in the model, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was conducted on the 22 variables. Following that, multivariate Cox proportional hazards regression analysis was applied to compute risk scores for β-alanine metabolism using the “predict” function in R. Subsequently, the samples were divided into high-risk and low-risk subgroups based on the median risk score of the train cohort. Survival analysis using the Kaplan-Meier method was then employed on the train, test1, test2, and test3 cohorts to assess the predictive capacity of the risk scores.

To ascertain the independent prognostic value of the risk score, both univariate and multivariate Cox regression analyses were executed. Only indicators that demonstrated statistical significance (p<0.05) in both analyses were considered as independent prognostic factors. Given the variations in clinical information between the TCGA and GEO databases, independent prognostic analyses were solely conducted in the TCGA cohort. To provide a quantitative approach for predicting the survival probability of HCC patients in clinical practice, we employed the “rms” package in R to construct a nomogram plot. This plot integrated independent prognostic indicators based on the TCGA cohort. Furthermore, we conducted a calibration plot to assess the alignment between actual and nomogram-estimated survival probabilities. Additionally, ROC curves were utilized to verify the diagnostic performance in predicting 1-, 3-, and 5-year survival rates.

Clinical significance and expression experimental verification of the EHHADH gene

The EHHADH gene constitutes a vital component of our HCC prognosis model. We have identified a significant correlation between this gene and various clinical-pathological indicators associated with HCC. As a result, our study exclusively delved into a thorough and comprehensive analysis, coupled with experimental validation of EHHADH. The findings presented in this study are predominantly derived from the BEST online platform, an open and encompassing multi-omics data visualization platform for diverse cancer types. Through queries involving HCC and EHHADH, we uncovered correlations between the gene and several key clinical indicators, including expression levels, age, grade, stage, BMI index, sorafenib sensitivity, microvascular invasion, among others. Furthermore, we conducted Kaplain-Meier survival analysis on multiple publicly available datasets to underscore the prognostic significance of EHHADH.

Utilizing the liver cancer tissue microarray (ZL-LivHcc962) procured from Shanghai Zhuoli Biotech Co., Ltd. (Zhuoli Biotech Co., Ltd., Shanghai, China), we conducted an analysis of EHHADH expression in both tumor tissues and their corresponding non-tumor counterparts. The tissue microarray underwent a 60-minute incubation in a constant temperature oven set at 59° C, followed by dewaxing and hydration using xylene and ethanol. Antigen retrieval was achieved through high-temperature and high-pressure treatment with EDTA. A 5% BSA blocking solution was administered dropwise and allowed to incubate at room temperature for 35 minutes. Subsequently, an optimal quantity of EHHADH antibody (Affinity Biosciences, DF4280, USA, 1:50) was introduced dropwise and incubated overnight at 4° C. The tissue microarray was then extracted, subjected to PBS washing, and exposed to a dropwise addition of the secondary antibody for a 35-minute incubation at 37° C. Diaminoaniline staining was employed for observation, followed by counterstaining with hematoxylin. The histochemistry score was computed based on the total area percentage and staining intensity, represented by the H-Score = ∑(pi×i) = (percentage of weak intensity×1) + (percentage of moderate intensity×2) + (percentage of strong intensity×3). Here, pi denotes the positive signal pixel area/cell count percentage, and i signifies the staining intensity. The H-score, ranging from 0 to 300, reflects a higher value corresponding to an augmented overall positive intensity of EHHADH.

Availability of data and materials

The datasets examined in this study are available in the Supplementary Materials or can be obtained by contacting the corresponding author.

Pan-cancer analysis of 22 βAMRGs

In order to explore the molecular characteristics of the 22 βAMRGs in diverse human cancers, comprehensive pan-cancer analyses were conducted, utilizing genomics, transcriptomics, and clinical data. Initially, we examined the mutational patterns of the 22 βAMRGs in human cancers and studied the prevalence of copy number variations (CNV). The results demonstrated that CNV occurred at a high frequency (ranging from approximately 10% to 80%) in various cancer types (Figure 1A). We also found the CNDP1 has a lower CNV frequency than other genes in human tumors. Then we researched the expression level of 22 βAMRGs in human tumors (Figure 1B). A significant number of the 22 βAMRGs exhibited differential expression between tumor tissues and adjacent normal samples in various human cancers. Furthermore, in nearly all human tumor tissues, the majority of βAMRGs were upregulated when compared to their expression levels in paired normal samples, particularly in cases of BRCA. We also found the SRM, GAD1 and SMS were significantly elevated in contrast with those in paired normal samples in human cancers. Next, we delved into the analysis of SNV within βAMRGs. Notably, BLCA, COAD, LUAD, LUSC, SKCM, STAD, and UCEC displayed a higher frequency of SNV, particularly in SKCM and UCEC (Figure 1C, 1D). Conversely, SNV frequency in other human tumors was relatively low, including ACC, CHOL, KICH, MESO, PCPG, TGCT, THYM, and UVM (Figure 1C, 1D). Remarkably, we observed that DPYD exhibited significantly higher SNV frequencies in human cancers, especially in SKCM, compared to other βAMRGs. These observations were visually illustrated through the heatmap and waterfall diagram of SNV (Figure 1C, 1D).

Figure 1. Pan-cancer analysis of 22 βAMRGs. (A) The frequency of CNV gains and losses was analyzed for the 22 βAMRGs across 20 different types of human cancers. The length of the line represents the variation frequency of the βAMRGs in pan-cancers. (B) The expression characteristics of the 22 βAMRGs were examined in 20 different types of human cancers, showing statistical significance (P < 0.05). (C, D) Heatmap and waterfall diagram were utilized to depict the SNV data of the 22 βAMRGs in pan-cancers. (E) Survival landscape analysis was conducted for the 22 βAMRGs across 20 different types of human cancers. Genes with P > 0.05 are represented by the white color, while red and blue colors indicate risk and protective genes, respectively. (F) The DNA methylation patterns of the 22 βAMRGs were examined in 20 different types of human cancers. The color gradient from orange to green signifies high to low methylation levels.

Subsequently, we explored the relationship of 22 βAMRGs expression with patient survival time (Figure 1E). By conducting a univariate Cox regression analysis, we identified βAMRGs that functioned as risk factors and those with a protective role (HR<1 and p<0.05). For instance, the most of the βAMRGs were the protective factors in the KIRC. It was worthy noticed that the ABAT βAMRG was the protective factors for human cancer. Methylation is a significant modification of proteins and nucleic acids, regulating gene expression and gene silencing. This epigenetic process plays a crucial role in various diseases, including cancer, aging, and Alzheimer’s disease, making it a vital focus of research in the field of epigenetics. DNA methylation is responsible for silencing specific genes, whereas demethylation leads to the reactivation and expression of these genes [20]. Hence, we researched that methylation patterns of βAMRGs in human cancer. According to the result (Figure 1F), the βAMRGs had complicated methylation patterns. For instance, the GAD2, AOC2 and DPYS had the hypermethylation in more than most of human cancers; on the contrary, most of βAMRGs showed the hypomethylation in more than most of human cancers. Interestingly, for LIHC, all the 22 βAMRGs showed the hypomethylation in human tumors.

Single-cell analysis of two HCC patients

According to the pre-established data quality control parameters, the single-cell quality control results for these two HCC patients were illustrated in Supplementary Figure 1A. The robust correlation between ncount_RNA and nfeature_RNA attested to the high quality of this batch of single-cell data (Supplementary Figure 1B). Notably, IGKC and IGHG1 genes emerged as significantly variable genes in this dataset (Supplementary Figure 1C). The single-cell atlas, both before and after batch correction using harmony, was presented in Supplementary Figure 1D, 1E. Overall, this study identified 31 cellular clusters (Figure 2A). Through automated annotation using singleR, seven distinct cell subgroups were conclusively recognized, encompassing T cells, B cells, monocytes, hepatocytes, smooth muscle cells, epithelial cells, and endothelial cells (Figure 2B, 2C).

Figure 2. Single-cell atlas of 2 HCC patients. (A) t-SNE dimensionality reduction subtypes; (B) t-SNE dimensionality reduction cell annotations; (C) singleR automated annotation results.

Regardless of the prediction algorithm employed, consistently obtained results indicate that beta-alanine metabolism was predominantly active in hepatocytes and epithelial cells (Figure 3A, 3B). The dimensionality reduction maps of beta-alanine metabolism features for each cell were depicted in Figure 3C, highlighting that the region housing hepatocyte clusters exhibited higher beta-alanine metabolism scores. To assess differences in beta-alanine metabolism between benign and malignant cells, the CopyKat algorithm was utilized for prediction (Supplementary Figure 2). As depicted in Figure 4A, the majority of hepatocytes were classified as malignant cells, while other cell types were identified as benign cells. This alignment with actual data suggests relatively accurate and credible predictions. Quantitative analysis results indicated elevated levels of beta-alanine metabolism in malignant cells (Figure 4B). The t-SNE dimensionality reduction map further elucidated that the region housing malignant cells exhibited heightened beta-alanine metabolism activity (Figure 4C). In conclusion, all findings consistently supported the study’s conclusions.

Figure 3. Assessing beta-alanine metabolism at the single-cell level. (A) Evaluation of the activation level of beta-alanine metabolism in each cell using six prediction algorithms; (B) Violin plots illustrating the variation in cellular activity levels of beta-alanine metabolism; (C) t-SNE plots depicting the dimensional distribution characteristics of beta-alanine metabolism.

Figure 4. Analysis of beta-alanine metabolism activity between benign and malignant cells. (A) t-SNE dimensional distribution characteristics of benign and malignant cells; (B) Differences in the activity of beta-alanine metabolism between benign and malignant cells; (C) t-SNE plots illustrating the distribution characteristics of beta-alanine metabolism between benign and malignant cells.

Cluster analysis for 343 HCC patients according to the βAMRGs scores

Taking into account the impacts of βAMRGs in HCC, we initially assessed the enrichment scores of βAMRGs using ssGSEA for a total of 343 patients. These HCC samples were then categorized into three distinct clusters (Cluster1, Cluster C2, and Cluster C3) based on the mRNA expression levels of βAMRGs (Figure 5A). Specifically, the clusters named C1, C2, and C3 comprised patients with inactive, normal, and active βAMRGs, respectively. The violin plot illustrated that C3 exhibited the highest enrichment score, while C1 had the lowest, followed by C2 (Figure 5B). Subsequently, we investigated the relationship between the three clusters and overall survival by plotting the survival curves. Among the clusters, C3 displayed the highest overall survival rates, whereas C1 exhibited the lowest (Figure 5C). Additionally, the overall survival rates for HCC patients in C2 were found to lie between those of C1 and C3, indicating that high βAMRGs scores acted as a protective indicator. The β-alanine metabolism played an essential role in the outcome of HCC, particularly in cases of high activity of β-alanine metabolism. In summary, based on these results, we conclude that the classification technique for HCC is accurate, reliable, and scientifically sound.

Figure 5. Cluster analysis for 343 HCC patients according to the βAMRGs scores. (A) The clustering of gene data reveals three distinct clusters, as depicted by the heatmap: Cluster 1(C1), Cluster 2 (C2), and Cluster 3 (C3), based on the levels of mRNA expression of βAMRGs. Cluster 2 (C2) comprises tumor patients with normal βAMRG expression, Cluster 3 (C3) includes those with active βAMRGs, while Cluster 1(C1) consists of individuals with inactive βAMRGs. (B) The violin plot illustrates the enrichment scores for the three clusters, arranged in the order of C3 > C2 > C1. The corresponding p-values are indicated above the respective clusters. (C) A survival curve represents the three distinct clusters. Cluster 3 exhibits the highest survival rate, while Cluster 1 demonstrates the lowest survival rate when compared among the three clusters. The x-coordinate represents survival time, while the y-coordinate represents survival rate.

Correlations of the βAMRGs scores with tumor related metabolic, immune and cell death pathways in the 3 clusters

An analysis using a heatmap was conducted to examine the activity of 42 conventional metabolic pathways within the three clusters of βAMRGs (Figure 6A). Significantly different activity levels were observed across most of the 42 pathways in the three clusters. Notably, the C3 cluster exhibited the highest activities of various metabolism pathways, including alanine aspartate and glutamate metabolism, beta alanine metabolism, and glycolysis/gluconeogenesis metabolism. Conversely, a few metabolic pathways such as ether lipid metabolism, riboflavin metabolism, and inositol phosphate metabolism exhibited the lowest activity in the C3 cluster. These findings aligned with our expectations, as β-alanine metabolism plays both positive and negative roles in tumor-related metabolic pathways. Additionally, the activity of 24 immune-related pathways across the three clusters was explored (Figure 6B). Interestingly, nearly all immune-associated pathways displayed greater activity in the C1 cluster compared to C2 and C3 clusters, including pathways such as B cell receptor signaling pathway, natural killer cell mediated cytotoxicity, and T cell receptor signaling pathway. This observation suggests that HCC malignancy is highest in the C1 cluster, or that patients within C1 exhibit immune escape, resulting in a poorer prognosis. Next, the activity of 10 cell death pathways was examined in the three clusters (Figure 6C). Significant differences were observed in nine out of the ten cell death pathways, with the exception of the PANoptosis cell death pathway. Notably, a positive correlation was observed between curroptosis and β-alanine activity.

Figure 6. Correlations of the scores of βAMRGs with tumor-related metabolic, immune, and cell death pathways in the three clusters were examined. (A) The activity levels of 42 conventional metabolic pathways were analyzed across the three clusters. (B) The activity levels of 24 immune-related pathways were assessed across the three clusters. (C) The activity levels of 10 cell death pathways were investigated across the three clusters. (* indicates p <0.05; ** indicates p < 0.01; *** indicates p < 0.001; **** indicates p < 0.0001).

Correlations between the βAMRGs scores and the immune cell infiltration in the 3 clusters

The tumor microenvironment (TME) encompasses a complex interplay among cells, extracellular matrix, and signaling molecules in the vicinity of tumor cells [21]. In the context of hepatocellular carcinoma (HCC), the TME assumes a critical role as an intrinsic factor contributing to the initiation, progression, invasion, and metastasis of the disease [22]. Particularly relevant to the pathogenesis of HCC, the TME serves as a pivotal regulator by delivering, inhibiting, or promoting growth signals. Consequently, it represents a valuable resource for identifying potential therapeutic targets, particularly pertaining to tumor-infiltrating immune cells (TIICs) [23]. And our comprehensive review of the existing literature did not reveal any prior investigations on the association between β-alanine and TIICs. In order to depict variations in immune cell infiltration levels among the C1, C2, and C3 subgroups, a heatmap was generated using TIMER, CIBERSORT, CIBERSORT−ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC algorithms (Figure 7A). The heatmap illustrating immune cell infiltration demonstrated that most immune cells exhibited the highest infiltration levels in the C1 cluster and the lowest levels in the C3 cluster. These results suggest that βAMRGs scores play a crucial role in modifying immune cell infiltration in HCC. Immune checkpoints are a molecular signaling pathway that regulates the immune response to maintain the balance of the immune system and prevent excessive activation. However, tumor cells can use these immune checkpoints to evade attack by the immune system, thereby promoting tumor growth and spread [24]. Hence, we studied the expression level of immune checkpoints related genes (ICGs) in different expression of βAMRGs clusters (Figure 7B). Our findings revealed a correlation between the βAMRGs-inactive cluster and the overexpression of immune checkpoint genes (ICGs), indicating that patients in the βAMRGs-inactive group exhibited the weakest anti-tumor immunity and poor prognostic outcomes.

Figure 7. The correlation between βAMRGs scores and immune cell infiltration within the three clusters was investigated. (A) The heatmap displays variations in immune cell infiltration levels across the three subgroups based on βAMRGs scores. (B) The box plot demonstrates discrepancies in immune checkpoint expression among the three subgroups. (C) The heatmap depicts the relationship between βAMRGs-associated genes and immune cell infiltration levels. (D) The plot shows the link between βAMRGs scores and immune cell infiltration levels. The size of the sphere on the right side represents the correlation strength (abs) and the color represents the corresponding p-value. (E–H) The scatter plot depicts the relationship of βAMRGs scores with 4 immune-infiltration-related substances. The βAMRGs scores were shown to have a positive link to the infiltration levels of Type I IFN Response and Type II IFN Response. The βAMRGs scores were shown to have a negative link to the infiltration levels of aDCs and Treg. (* indicates p <0.05; ** indicates p < 0.01; *** indicates p < 0.001; **** indicates p < 0.0001).

Our study also endeavors to elucidate the relationship between βAMRGs’s levels and immune microenvironment (Figure 7C). We observed a predominantly negative correlation between the majority of βAMRGs (e.g., UPB1, EHHADH, ALDH9A1, ALDH7A1) and immunocyte infiltration, as well as immune-related functions. Conversely, βAMRGs such as SRM, SMS, GAD1, DPYD, and AOC3 exhibited a positive association with the levels of immune cell infiltration. And with the exception of SRM, SMS, GAD2, and GAD1, all βAMRGs exerted a positive influence on the Type II IFN Response. Finally, we evaluated the relationship between βAMRGs scores and immune cell infiltration (Figure 7D–7H). Notably, mast cells, neutrophils, NK cells, Type I IFN Response, and Type II IFN Response demonstrated a positive correlation with βAMRGs scores, while the remaining infiltrating immune cells exhibited a negative correlation, especially for aDCs and Treg.

Construction and verification of a novel βAMRGs related prognostic signature(βAMRGs-RPS) for predicting the clinical outcomes of patients with HCC

A total of 22 βAMRGs from previous studies were utilized for the LASSO-Cox regression analysis. Subsequently, seven out of the 22 βAMRGs were selected for further multivariate Cox regression analysis (Supplementary Figure 3A, 3B). Finally, a multivariate Cox proportional hazards regression analysis was conducted, integrating six central hub βAMRGs (i.e., SMS, HIBCH, GAD1, GAD2, CNDP1, EHHADH), to construct the βAMRGs-RPS (Supplementary Figure 3C). In cohort 1, which comprised 343 TCGA samples, the three GEO cohorts were merged to form cohort 2. Cohort 1 samples were randomly divided in a 6:4 ratio, with 60% assigned to the training group, 40% to the test1 group, and all samples of TCGA assigned to the test2 group. The cohort 2 samples constituted the test3 group. Based on the median risk score, the patients in the training set were categorized into high- and low-risk groups. The expression levels of the six genes were visualized using a heatmap, depicting the differences between the high- and low-risk score groups (Figure 8A). The survival curve analysis revealed that the high-risk group had a poorer overall survival compared to the low-risk group (Figure 8B). Furthermore, the HCC patients were stratified into low- and high-risk subgroups based on the median risk score (Figure 8C). The high-risk subgroup exhibited a significantly higher mortality rate compared to the low-risk subgroup, as indicated by the risk score distributions and survival status (Figure 8D). These findings indicate that the βAMRGs-RPS can accurately distinguish HCC patients with a favorable or unfavorable prognosis.

Figure 8. Establishment of a novel βAMRGs-RPS in the train cohort. (A) Heatmap illustrating the expression levels of six genes between subgroups with high and low risk scores. (B) Prognostic prediction depicted by the Kaplan-Meier survival curve comparing subgroups with high and low risk scores. (C) Risk score curve plot exhibiting the distribution of individual risk scores, with the patients categorized into low-risk (green) and high-risk (red) groups. (D) Risk score scatter plot, where red dots indicate deceased patients and green dots indicate surviving patients.

To validate the reliability and stability of the βAMRGs-RPS, heatmaps were generated for the test1, test2, and test3 cohorts, demonstrating a similar expression pattern of the six βAMRGs as observed in the training cohort (Supplementary Figures 4A, 5A, 6A). Notably, both internal and external validation of the survival data demonstrated that HCC patients with low-risk scores exhibited significantly higher overall survival rates (p < 0.05) (Supplementary Figures 4B, 5B, 6B). Additionally, the HCC samples in the test1, test2, and test3 cohorts were also classified into low- and high-risk populations using the same bioinformatics methodology (Supplementary Figures 4C, 5C, 6C). It is important to note that the median risk score from the training cohort served as a consistent benchmark for separating HCC samples. The risk score distributions and survival status in the internal validation cohorts (test1 and test2) as well as the external validation cohort (test3) showed similar trends to those observed in the training cohort (Supplementary Figure 4D, 5D, 6D). These results highlight the reliability and accuracy of the βAMRGs-RPS.

Constructing a nomogram plot by integrating independent prognostic indicators based on the TCGA cohort

To further investigate the precision of this signature, a comprehensive analysis was undertaken to evaluate its potential as an independent prognostic factor for overall survival (OS) within the TCGA cohort. Both univariate and multivariate analyses were conducted, incorporating the risk score and clinical traits as variables. The results of the univariable Cox regression analysis revealed significant correlations between OS and cancer status, stage, and the risk score (p < 0.05) (Figure 9A). Subsequently, the multivariate analysis confirmed cancer status, stage, and the risk score as independent predictors for OS (p < 0.05) (Figure 9B). Based on these three variables, a nomogram for OS was constructed (Figure 9C). The calibration curves demonstrated a high degree of concordance between the actual proportions of 1-, 3-, and 5-year OS and the predicted probabilities derived from the nomogram (Figure 9D). Furthermore, to assess the accuracy and stability of the nomogram, a receiver operating characteristic (ROC) curve analysis was performed. The AUC values for the 1-, 3-, and 5-year ROC curves were found to be 0.740, 0.718, and 0.751, respectively (Figure 9E). The ROC curve further substantiates the novelty, reliability, and practicality of our findings.

Figure 9. Establishment of nomogram. (A) Univariate analyses were conducted to assess the relationship between risk scores and relevant clinical parameters and OS in the TCGA cohort. (B) Multivariate analyses were performed to evaluate the relationship between risk scores, relevant clinical parameters, and OS in the TCGA cohort. (C) Development of a signature-based prognostic nomogram was undertaken to predict OS in HCC. (D) Calibration curves of the nomogram prediction of 1-year, 3-year, and 5-year OS of patients in TCGA cohort. (E) ROC curve of the risk score. The AUC values for the 1-, 3-, and 5-year ROC curves were found to be 0.740, 0.718, and 0.751, respectively.

Clinical significance and expression experimental verification of the EHHADH gene

The EHHADH gene played a pivotal role in our prognostic model, and we have identified it as a potential protective protein in HCC. Our findings reveal a significant decrease in EHHADH gene expression in tumor tissues compared to adjacent non-cancerous tissues (Figure 10A–10C). Moreover, patients with a younger age and lower BMI index exhibit reduced levels of EHHADH gene expression (Figure 10D, 10E). As histological and pathological stages progress, the expression of the EHHADH gene further diminishes, strongly indicating its potential as a protective factor in HCC (Figure 10F–10I). Additionally, patients with microvascular invasion and those with a poor response to sorafenib treatment also exhibit lower levels of EHHADH gene expression (Figure 10F–10I). Decreased levels of EHHADH gene expression often correlate with a poorer prognosis for HCC patients (Figure 10J–10O), and targeting the activation of EHHADH may represent a viable strategy to prolong patients’ survival. Finally, we complemented the immunohistochemical experiments of EHHADH with a tissue microarray of liver cancer (Figure 11). The results illustrated a clear trend of downregulation in the protein levels of EHHADH in the majority of liver cancer samples.

Figure 10. Association of EHHADH expression with clinical traits. The expression traits of EHHADH in (A) GSE144269, (B) GSE54236, (C) TCGA cohorts. Association of EHHADH expression with (D) age, (E) BMI, (F) grade, (G) microvascular invasion, (H) stage, and (I) sorafenib response. Prognostic performances of EHHADH in (J–O) multiple HCC cohorts.

Figure 11. Immunohistochemical experiment and statistical analysis of the EHHADH gene based on HCC tissue microarray. (A–H) merely represent the rows of liver cancer tissue chips, without any special significance that needs to be reflected in the paper.