Comprehensive molecular characteristics of hepatocellular carcinoma based on multi-omics analysis

Abstract

Background

The high heterogeneity of hepatocellular carcinoma (HCC) poses challenges for precision treatment strategies. This study aims to use multi-omics methodologies to better understand its pathogenesis and discover biomarkers.

Methods

Quantitative proteomics was used to investigate hepatocellular carcinoma tissues (HCT) and their corresponding adjacent non-tumor tissues (DNT), obtained from six HCC patients. Untargeted metabolomics was applied to analyze the metabolic profiles of HCT and DNT of ten HCC patients. Statistical analyses, such as the Student’s t-test, were performed to identify differentially expressed proteins (DEPs) and metabolites (DEMs) between the two groups. The functions and metabolic pathways involving DEPs and DEMs were annotated and enriched using the gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) databases. Bioinformatics methods were then utilized to analyze consistency between proteomics and metabolomics results, leading to identification of potential biomarkers along with key altered pathways associated with HCC.

Results

This study identified 1556 DEPs between HCT and DNT samples. These DEPs were primarily enriched in crucial biological pathways such as amino acid degradation, fatty acid metabolism, and DNA replication. Subsequently, the analysis of metabolomics identified 500 DEMs that mainly participated in glycerophospholipid metabolism, the phospholipase D signaling pathway, and choline metabolism related to cancer. Integrated analysis of proteomics and metabolomics data unveiled significant dysfunctions in bile secretion, multiple amino acid and fatty acid metabolic pathways among HCC patients. Further investigation revealed that five proteins (PTP4A3, B4GALT5, GAB1, ME2, and PKM) along with seven metabolites (PI(6 keto-PGF1alpha/16:0), 13, 16, 19-docosatrienoic acid, PA(18:2(9Z, 12Z)/20:1(11Z)), Citric Acid, PG(20:3(6, 8, 11)-OH(5)/18:2(9Z, 12Z)), Spermidine, and N2-Acetylornithine) exhibited excellent diagnostic efficiency for HCC and could serve as its potential biomarkers.

Conclusion

Our integrated proteome and metabolome analysis revealed 10 key HCC-related pathways and proposed 12 potential biomarkers, which may enhance our understanding of HCC pathophysiology and be helpful in facilitating early diagnosis and treatment strategies.

Introduction

Hepatocellular carcinoma (HCC) is the predominant histological subtype of primary liver cancer, ranking sixth in global malignant tumor incidence and third as a leading cause of cancer-related mortality [1]. Annually, there are approximately 740,000 new cases of HCC and 640,000 deaths worldwide [2]. Despite significant advancements in the diagnosis, treatment, and prevention of HCC in recent years, its intricate pathogenesis and diverse clinical manifestations remain focal points and challenges for scientific research. The majority of HCC patients are diagnosed at an advanced stage where radical treatment options are not feasible. Studies have revealed that the 5-year survival rate of HCC patients undergoing surgical resection in the intermediate and advanced stages ranges from only 50–68%, with a recurrence rate as high as 50-70%. Conversely, HCC patients diagnosed and treated at an early stage exhibit a more than 70% 5-year survival rate, accompanied by a mere 10% recurrence rate [3, 4]. Consequently, employing advanced research methodologies for an in-depth exploration of HCC pathogenesis and the identification of novel diagnostic markers holds immense significance in enhancing early detection rates, overall survival rates, and quality of life among HCC patients.

In recent years, the rapid development of post-genomic technology has made omics research methods an indispensable tool for studying disease pathogenesis due to their multidimensional, systematic, and diverse advantages. Proteomics enables comprehensive analysis of protein expression, structure, and function in organisms to reveal dynamic changes and regulatory networks [5]. Metabolomics provides a deep understanding of cellular metabolism changes and their relationship with disease occurrence and development by studying all metabolic substances and pathways in an organism [6]. Using proteomics methods to analyze liver tissue and plasma of HCC patients, researchers identified alterations in spliceosome pathway activation, cholesterol homeostasis disorder, and up-regulation of coagulation and complement related pathways in HCC patients. Additionally, it is suggested that aldo-keto reductase family 1 member B10 (AKR1B10), cathepsin A (CTSA), sterol o-acyltransferase1 (SOAT1), pyrroline-5-carboxylate reductase 2 (PYCR2), and alcohol dehydrogenase 1 A (class I), alpha polypeptide (ADH1A) have potential as biomarkers for the diagnosis and prognosis of HCC [ 7,8,9]. Based on metabolomics methods, researchers investigated the plasma and liver tissues of HCC patients and observed abnormalities in metabolic pathways including tricarboxylic acid cycle (TCA cycle), glycolysis, purine metabolism, fatty acid metabolism, and retinol metabolism. Free fatty acids such as betaine and propionylcarnitine, retinol derivatives like retinaldehyde, uric acid levels along with bile acids and elaidic acid were found to play significant roles in these metabolic pathways. These metabolites hold promise as potential diagnostic markers for HCC [10,11,12,13,14,15,16]. The application of liver tissue proteomics and metabolomics is particularly crucial in the study of HCC. Utilizing paired omics analysis of different tissues in patients not only effectively mitigates individual differences such as age, gender, and region but also significantly enhances the accuracy of screening organ-specific biomarkers. This approach enables a more accurate reflection of the physiological and pathological state of the liver itself, independent from other organs or systemic influences [6].

With the advancement of science and technology, multi-omics integrated analysis has emerged as a more precise, comprehensive, and efficient approach for disease research. By integrating transcriptomics, proteomics, and metabolomics analyses, researchers have conducted thorough investigations into the pathogenesis of malignant tumors such as liver cancer, colorectal cancer, and breast cancer. Consequently, they have identified potential diagnostic markers, novel therapeutic targets, and prognostic indicators. This innovative methodology offers enhanced accuracy and effectiveness in the prevention, diagnosis, treatment, and prognosis of these diseases [17,18,19]. Despite the progress made, there remains a dearth of comprehensive exploration into the pathogenesis and biomarkers of HCC through the integration of liver tissue proteomics and metabolomics. Moreover, the overall research status concerning protein and metabolite regulatory networks in HCC patients remains unclear, with highly specific diagnostic markers yet to be identified. Therefore, this study aims to comprehensively apply liver tissue proteomics and metabolomics technologies to compare and analyze paracancerous liver tissues with HCC tissues’ protein and metabolite expression profiles. The objective is to identify differentially expressed proteins and metabolites closely associated with HCC while investigating their correlation and role in its pathogenesis. Ultimately, this study aims to provide novel insights for early diagnosis and personalized treatment of HCC.

Materials and methods

Participants and sample collection

Ten patients with HCC who were admitted to the Second Hospital of Lanzhou University between March 2023 and April 2024 were included in this study. Non-targeted metabolomics analysis was conducted on selected samples of hepatocellular carcinoma tissues (HCT) and their corresponding adjacent non-tumor tissues (DNT). All patients had a confirmed pathological diagnosis of HCC and had not undergone any prior treatment, including chemotherapy, radiotherapy, immunotherapy, or other therapies. The metabolic profile changes in ten HCC patients were analyzed using the non-targeted metabolomics approach. Additionally, the data-independent acquisition (DIA) quantitative proteomics method was employed to further investigate the protein and metabolomic profiles of six HCC patients by analyzing their HCT and DNT samples jointly. The HCT samples were collected from the central region of surgically resected liver tissue containing solid tumors, while DNT samples were obtained from an area located 5 cm outside the boundary of the solid tumor. Immediately after surgical resection, all tissue specimens were flash-frozen in liquid nitrogen to preserve their biological activity and stored at -80 °C until analysis.

DIA quantitative proteomic analysis

Extraction of protein

Under frozen conditions, excise tissue weighing 100 ± 0.5 mg and transfer it to a grinding tube. Add an appropriate volume of protein lysis buffer (8 M urea + 1% SDS) supplemented with protease inhibitors. Agitate the mixture three times using a cryogenic grinder for 40 s each time. Place the mixture on ice for 30 min, vortexing for 5–10 s every 5 min to facilitate lysis. Subsequently, centrifuge at 4 °C, 16,000 g for 30 min and collect the supernatant. Quantify the protein content in the supernatant using Bradford assay and sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE).

Peptide sample preparation and DIA mass spectrometry detection

Firstly, a protein sample of 100 µg was taken and mixed with lysis buffer and triethylammonium bicarbonate buffer (TEAB) to achieve a final concentration of 100 mM. Tris (2-carboxyethyl) phosphine (TCEP) was then added to reach a final concentration of 10mM, and the reaction was carried out at 37 °C for 60 min. Subsequently, Iodoacetamide (IAM) was added to attain a final concentration of 40mM, and the mixture was incubated in darkness at room temperature for 40 min. Afterwards, the mixture was subjected to centrifugation after adding six volumes of pre-cooled acetone and allowing it to stand at -20 °C for four hours in order to separate the precipitate. The precipitate was dissolved in 100µL of TEAB with a concentration of 100mM and digested overnight at 37 °C using trypsin at a mass ratio of 1:50. The resulting peptide mixture was then dried under vacuum to remove solvents, followed by resuspension in trifluoroacetic acid (TFA) solution with a concentration of 0.1%. Desalting using an HLB column was performed to eliminate impurities from the desalted peptides which were subsequently dried again under vacuum to remove any remaining solvents. Finally, quantification of peptides took place using NANODROP ONE UV-Visspectrophotometer (Thermo Fisher Scientific, San Jose, USA).

Equal amounts of peptides were initially dissolved in a loading buffer for mass spectrometry and subsequently injected into a vanquish Neo chromatography system (Thermo Fisher Scientific, San Jose, USA). High-performance liquid chromatography (HPLC) separation was conducted using a uPAC High Throughput column (75 μm × 5.5 cm, Thermo Fisher Scientific, San Jose, USA). A gradient elution method was employed, utilizing mobile phase A containing 2% acetonitrile and 0.1% formic acid, along with mobile phase B consisting of 80% acetonitrile and 0.1% formic acid. During an 8-minute chromatographic separation, the ratio between mobile phases A and B varied according to a predefined gradient to optimize sample separation. Subsequent to the separation process, the peptides entered an astral mass spectrometer (Thermo Fisher Scientific, San Jose, USA) for ionization and measurement of their mass-to-charge ratios, thereby generating mass spectrometry data. Data acquisition and initial processing were managed by Thermo xcalibur software (Thermo Fisher Scientific, San Jose, USA), employing a SWATH-based DIA approach throughout the entire procedure.

Differentially expressed proteins analysis

The proteome discoverer software (Thermo Fisher Scientific, San Jose, USA) was utilized for hierarchical database construction and data analysis. Subsequently, the database search results were imported into Spectronaut™ software (Biognosys AG, Schlieren, Switzerland) to establish the spectral library. The Spectronaut™ software was employed to extract DIA raw data peaks based on ion pair information from the spectrum library. Default settings of the software and iRT correction were applied during analysis. Up to 6 specific peptides per protein and up to 6 child ions per peptide were selected. Protein quantification results underwent preprocessing by eliminating duplicate values and filling missing values. Differential expression proteins (DAPs) were identified using univariate analysis and other methods with a screening threshold of the absolute value of the logarithm base 2 of the fold change (|Log₂FC|) > = 1 and adjusted p-value < 0.05.

Metabolomics analysis

Metabolite extraction

The tissue sample weighing 100 mg was placed in a 2mL centrifuge tube, followed by the addition of a ground bead with a diameter of 6 mm. Subsequently, 400µL of an extraction solution consisting of methanol and water in a volume ratio of 4:1, containing an internal standard L-2-chlorophenylalanine at a concentration of 0.02 mg/mL, was added to the centrifuge tube for metabolite extraction. Next, the sample solution was homogenized using a frozen tissue grinder for 6 min (-10 °C, 50 Hz) and then subjected to low temperature ultrasound extraction for 30 min (5 °C, 40 kHz). After sonication, the centrifuge tube was cooled at -20 °C for 30 min and subsequently centrifuged for another 15 min (4 °C, 13000 g). Finally, the supernatant was transferred into an injection vial equipped with an inner cannula for analysis.

LC-MS/MS analysis

The LC-MS analysis was performed using an ultra-high performance liquid chromatography system coupled with a Q Exactive hybrid mass spectrometer (UHPLC-Q Exactive). Chromatographic conditions: 2µL of the sample was injected onto an HSST3 column for separation prior to mass spectrometric detection. Mobile phase A consisted of 95% water and 5% acetonitrile (containing 0.1% formic acid), while mobile phase B was a mixture of 47.5% acetonitrile, 47.5% isopropanol, and 5% water (also containing 0.1% formic acid). The separation gradient started with mobile phase B linearly increasing from 0 to 5% in 0–0.1 min, then from 5 to 25% in 0.1–2 min. From 2 to 9 min, B linearly rose to 100%, staying constant at 100% till 13 min. In 13.0-13.1 min, B linearly decreased to 0%. The flow rate was set at 0.40 mL/min, and the column temperature was maintained at 40 °C during the experiment. The mass spectrometric data was collected using a Thermo UHPLC-Q Exactive mass spectrometer equipped with an electrospray ionization (ESI) source operating in either positive or negative ion mode. The optimal conditions were set as follows: heater temperature, 400 °C; capillary temperature, 320 °C; sheath gas flow rate, 40arb; auxiliary gas flow rate, 10arb; ion-spray voltage floating (ISVF), -2800 V in negative mode and 3500 V in positive mode, respectively. Data acquisition was performed in Data-Dependent Acquisition (DDA) mode with the following settings: normalized collision energy, 20-40-60 V rolling for MS2; isolation window, 1.2 m/z; maximum injection time, 100ms; AGC target, 3e6; dynamic exclusion, 10s. The full MS resolution was set to 70,000, and the MS2 resolution was set to 17,500. The detection was carried out over a mass range of 70–1050 m/z for both MS1 and MS2 scans.

Differentially expressed metabolites analysis

The raw LC-MS data were imported into Progenesis QI (Waters Corporation, Milford, USA) for extensive preprocessing, including baseline filtering, peak identification, integration, retention time correction, and peak alignment. This comprehensive workflow resulted in a structured data matrix, consisting of retention times, mass-to-charge ratios (m/z), and peak intensities. Subsequently, the MS spectra within this data matrix were matched against both public database and Majorbio’s proprietary database to identify and annotate the metabolites presented in the samples. After searching the database, the matrix data was preprocessed, and the 80% rule was used to remove the missing values [20]. In order to reduce the error caused by sample preparation and instrument instability, the sum normalization method was used to normalize the response intensity of the spectral peak of the sample, and the normalized data matrix was obtained. At the same time, the variables with relative standard deviation (RSD) > 30% of quality control samples were deleted, and log10 logarithmic processing was performed to obtain the data matrix for subsequent analysis [20, 21]. Principal component analysis (PCA), partial least squares discriminant analysis (PLS-DA), and univariate analysis were utilized to screen for differentially expressed metabolites (DEMs) between groups. DEMs were identified based on the following criteria: FC value ≥ 1, variable importance in projection (VIP) score from the PLS-DA model ≥ 1, and a statistical significance level with adjusted p-value < 0.05 [22, 23].

Statistical analysis

Statistical analyses were performed using R (version 4.4.1). Normality was assessed using the Shapiro-Wilk test from the “stats” package, and group comparisons were conducted using Student’s t-test. For multivariate statistical analysis, PCA and PLS-DA were performed using the “ropls” packages (version 1.6.2). Data were scaled (unit variance scaling and mean centering) for PCA, and PLS-DA components were determined by 10-fold cross-validation. Spearman correlation analysis was conducted using the “stats” package. Data visualization was performed with “ggplot2” package (version 3.4.0). The gene ontology (GO) and human metabolome database (HMDB) databases were used for functional annotation analysis of all DEPs and DEMs. Kyoto encyclopedia of genes and genomes (KEGG) database metabolic pathways involved in the differences between proteins and metabolites annotation analysis. GO enrichment analysis was conducted with Goatools (version 1.4.4), using the Fisher exact test with Benjamini-Hochberg false discovery rate (FDR) correction (adjusted p-value < 0.05 considered significant). KEGG pathway enrichment analysis was performed using the “scipy” package (version 1.7.3) in Python (version 3.8), with the Fisher exact test and FDR correction applied. Relevant networks were constructed And visualized using Cytoscape software (version 3.7.2). All analyses were performed on the Majorbio cloud platform (cloud.majorbio.com). All statistical tests were two-sided, and adjusted p-value < 0.05 was considered statistically significant.

Result

Proteomics

Screening of deps

In this study, DIA quantitative proteomic analysis was performed on HCT and DNT tissues from six HCC patients, and a total of 126,081 peptides were identified and 8435 proteins were quantified. Based on the screening criteria (|Log₂FC| >=1 and adjusted p-value < 0.05), a total of 1556 DEPs were identified, of which 1148 proteins showed an up-regulated trend and 408 proteins showed a down-regulated trend in HCC patients (Fig. 1A). The results of cluster analysis showed that there were significant differences in DEPs distribution between HCT and DNT, and the intra-group reproducibility was relatively stable (Fig. 1B).

Fig. 1

Characterization of differentially expressed proteins. (A) Volcano plot illustrating the differential expression of proteins. (B) Hierarchical clustering heatmap depicting the differential protein profiles. The abbreviation NC denotes paracancerous tissue, while Ca refers to hepatocellular carcinoma tissue

Functional analysis of deps

In order to deeply understand and analyze the specific functions and roles of DEPs in biological processes, we performed detailed subcellular localization prediction for 1556 DEPs, and performed comprehensive protein function annotation and enrichment analysis based on GO and KEGG databases. The results of subcellular localization showed that the DEPs were mainly concentrated in the cytoplasm (41%), nucleus (20%) and mitochondria (15%) (Fig. 2A). The GO annotation results showed that the top 20 GO terms were mainly focused on biological process (BP) such as cellular process, metabolic process and biological regulation, followed by molecular function (MF) such as protein binding, catalytic activity and molecular function regulation activity. While the cellular component (CC) contains cellular anatomical entities and protein complexes (Fig. 2C). KEGG annotation results showed that 50% of the top 20 KEGG metabolic pathways in abundance belonged to the metabolic category, and among these metabolic pathways, the number of DEPs involved in amino acid metabolism, carbohydrate metabolism, and lipid metabolism was particularly significant (Fig. 2D). According to the adjusted p-value, GO and KEGG enrichment analysis were performed on the expression of DEPs. GO enrichment analysis showed that the DEPs were mainly enriched in cytochrome P450 epoxidase pathway, primary alcohol metabolism process and glutamine family amino acid metabolism process. In terms of cellular components, it mainly involves DNA prereplication initiation complex, CMG complex, mitochondrial large ribosomal subunit and organelle large ribosomal subunit. In terms of molecular functions, the main activities were arachidonic acid monooxygenase activity, arachidonic acid cyclooxygenase activity, aromatase activity and oxidoreductase activity (Fig. 2E). KEGG enrichment analysis showed that the DEPs were mainly enriched in valine, leucine and isoleucine degradation, fatty acid degradation, DNA replication and retinol metabolism pathways (Fig. 2F).

In order to determine the interaction between DEPs, we used the STRING (https://string-db.org/) database to predict the functional correlation between two groups of DEPs from the genomic correlation of their coding genes. The top 25 proteins with the highest connectivity degree were selected to draw the protein-protein interaction (PPI) network diagram (Fig. 2B). All the 25 proteins showed a trend of up-regulation in HCC, and the 10 proteins with the lowest adjusted p-value and high connectivity were POLR1B, MRTO4, GTPBP4, NIFK, TYMS, POLR1A, MRPL4, CDK1, BRIX1, and GRWD1 (Table 1).

Fig. 2

Functional analysis of differentially expressed proteins. (A) Subcellular localization mapping of differentially expressed proteins. Abbreviations: PM, plasma membrane; ER, endoplasmic reticulum; EX, extracellular; GA, golgi apparatus; PE, Periplasm; LY, lysosomal. (B) PPI network analysis of differentially expressed proteins. The intensity of the color of the circle in the figure corresponds to the number of nodes connected by the protein; a darker shade indicates a higher number of nodes connected by the protein. (C) GO term analysis for the differentially expressed proteins in HCC. (D) KEGG annotation for the differentially expressed proteins in HCC. (E) GO enrichment analysis for the identification of enriched terms among the differentially expressed proteins. The abscissa represents the enrichment rate, which signifies the ratio of enriched proteins in a specific GO term to the total number of proteins annotated to that term. A higher ratio indicates a greater level of enrichment. The ordinate denotes the GO term itself. Each bubble in the figure corresponds to a distinct GO Term, and its size is proportional to the number of concentrated proteins enriched by that particular term. The varying colors of the bubbles represent adjusted p-values, with smaller values indicating redder hues. Abbreviations: CC, cellular component; BP, biological process; MF, molecular function. (F) KEGG pathway enrichment analysis to identify significantly enriched pathways among the differentially expressed proteins. Abbreviations: M, metabolism; GIP, genetic information processing; HD, human diseases

Table 1 The PPI analysis identified the top 25 proteins with the highest degree of connectivity

Metabolomics

Screening of DEMs

In this study, PCA and PLS-DA analysis were used for pattern recognition of HCT and DNT in 10 HCC patients, respectively. PCA and PLS-DA analysis showed that there was a good separation trend between HCT and DNT (Fig. 3A and B). In order to judge the model quality, 200 permutation tests were performed on the PLS-DA model, and the results showed that the model was good and did not overfit (Fig. 3C). According to the screening conditions (FC ≥ 1, VIP of the PLS-DA model ≥ 1, and adjusted p-value < 0.05), a total of 500 DEMs were screened in this study, of which 195 showed an up-regulated trend and 305 showed a down-regulated trend in HCC patients (Fig. 3D). The HMDB compound classification results showed that the DEMs were classified into 11 categories. Among them, lipids and lipoid molecules (189), organic acids and their derivatives (95), organic heterocyclic compounds (66), and organic oxygen-containing compounds (35) were in the majority (Fig. 4).

Fig. 3

Characterization of metabolomics between HCT and DNT. (A) Plot of PCA scores between the two groups. After conducting dimension analysis, the reduced samples were plotted as relative coordinate points on the principal components p1 and p2. The distance between each point reflects the level of aggregation and dispersion among samples, with shorter distances indicating higher similarity and longer distances suggesting greater dissimilarity. (B) Plot of PLS-DA scores between the two groups. The PLS-DA score plot is commonly used to visually demonstrate the model’s classification effectiveness, where a larger separation between the two sample groups indicates a more pronounced classification effect. (C) PLS-DA permutation test plot. The abscissa represents the permutation retention of the test, while the ordinate represents R² (red dot) and Q² (blue triangle) values. The dashed lines correspond to regression lines for R² and Q² respectively. Notably, the regression line for R² is positioned above that of Q². Moreover, the intercept between the regression line for Q² and y-axis is -0.0926, indicating a well-fitted model with high predictability suitable for subsequent data analysis. (D) Volcano plot of differentially expressed metabolites. The horizontal axis represents the log₂FC in expression between the two groups, while the vertical axis represents the -log10(p-value) of statistical test for metabolite expression difference. Higher values indicate greater significance in expression differences. Each dot on the plot corresponds to a specific metabolite, with dot size indicating VIP value. Red dots represent up-regulated metabolites, blue dots represent down-regulated metabolites, and gray dots represent non-differential metabolites. The position of each dot reflects its level of significance in terms of expression difference

Fig. 4

The classification of compounds in HMDB

Among the top 20 DEMs with the lowest adjusted p-value, 10 were up-regulated in HCC (PI (6 keto-PGF1alpha/16:0), 13,16, 19-docosatrienoic acid and PA (18:2(9Z,12Z)/20:1(11Z)), etc.), and 10 were down-regulated (N-Amidino-L-aspartate, Lupex and 2,3,4, 5-tetrahydro-2-pyridinecarboxylic acid, among others) (Table 2). At the same time, correlation heat map showed that most of the above 20 DEMs had obvious positive correlation (Fig. 5A). PLS-DA was used as the supervised model, and 7-fold cross validation was used to test the prediction of different changes in pairing. The VIP score analysis of the first principal component identified the top 20 metabolites differentially associated with HCC occurrence. Among these metabolites, 4 were found to be up-regulated in HCC: D-glucurono-6,3-lactone, Pectic acid, Tryptophyl-Valine, and Senkyunolide N. Additionally, 16 of them were down-regulated, including galactosyl lactate, 2,3,4,5-tetrahydro-2-pyridinecarboxylic acid, and N, N-dimethyl-N’-(8-methyl-5 H-pyrimido[5,4-b] indol-4-yl) ethane-1,2-diamine (Fig. 5B).

Table 2 The top 20 metabolites with significant differences between HCT and DNT

Fig. 5

The metabolite features with the top 20 differences between HCT and DNT. (A) A correlation heatmap was generated for the top 20 differentially expressed metabolites. (B) VIP value analysis was performed on the top 20 differential metabolites. On the left side, a dendrogram of metabolite clustering is presented, where closer branches indicate similar expression patterns among all metabolites in the samples. Each column represents a sample. Each row represents a metabolite, and color indicates its relative expression level in this group of samples. The gradient color block provides a visual representation of the correspondence between color gradient and numerical size. On the right side, a VIP bar chart displays the contribution value of each metabolite to the difference between two groups (threshold set at 1). A longer bar indicates a greater difference in abundance between HCT and DNT groups. The color of each bar reflects significance levels (adjusted p-value) for differences between groups: darker colors represent smaller adjusted p-values (-log10(adjusted p-value)). Additionally, asterisks (*) denote P < 0.05, ** denotes P < 0.01, and *** denotes P < 0.001

Discussion

The objective of this study is to investigate the potential pathogenesis and biomarkers of HCC by conducting a combined omic analysis on HCT and DNT tissues from HCC patients. By performing functional and association analyses on DEPs and DEMs, along with their biological significance, this study identified 5 proteins such as PTP4A3 and 7 metabolites such as PI (6 keto-PGF1alpha/16:0) that may assist in the diagnosis of HCC. Furthermore, metabolic pathways such as tryptophan metabolism, butanoate metabolism, pyruvate metabolism, and abnormal bile secretion were found to be closely associated with symptoms like immune escape, abnormal energy metabolism, and lipid peroxidation in HCC patients.

The proteins PTP4A3, B4GALT5, GAB1, ME2, and PKM were found to exhibit high expression levels in HCC and demonstrated good diagnostic potential for HCC (AUC > 0.90). PTP4A3 is a phosphatase that plays a crucial role in cellular signaling by catalyzing the dephosphorylation of tyrosine or serine/threonine residues in protein molecules [24]. Its involvement in the development, metastasis, and recurrence of various cancers has been confirmed, particularly in solid tumors such as colorectal cancer, breast cancer, and HCC where its overexpression often indicates a poor prognosis [25, 26]. In our study, PTP4A3 showed a significant difference between HCT and DNT, with an adjusted p-value < 0.01 Correlation analysis further revealed a strong positive correlation between PTP4A3 and key metabolites 13-16-19-docosatrienoic acid and citric acid. These findings collectively emphasize the important role of PTP4A3 in HCC development and progression. The combined analysis of proteomics and metabolomics in this study revealed that B4GALT5 exhibited the highest FC (FC: 926.6) among the five potential protein biomarkers identified in this study. It was also the only one found to be significantly associated with PG (18:2 (z, z) / 18:1 (z) 9), PG (David,8,11(6)-OH no more(5)/(z, z)), PI (6 keto-PGF1alpha/16:0), and PA (18:2(9Z,12Z)/20:1(11Z)). These findings were consistent with Luo et al.‘s research [27], which demonstrated that B4GALT5 was highly expressed in HCC patients and enriched in lipid metabolism-related pathways. Additionally, KEGG enrichment analysis indicated that B4GALT5 primarily participated in sphingolipid metabolism and mucin type O-glycan biosynthesis pathways. This suggests a potential role for B4GALT5 in promoting HCC cell proliferation, metastasis, invasion, and ultimately contributing to poor prognosis among HCC patients [28]. As a versatile adaptor protein, GAB1 can respond to various stimuli such as extracellular growth factors, cytokines, and immune cell receptors. It is activated by multiple kinases, thereby regulating tumor proliferation, cell cycle progression, evasion of apoptosis, invasion and migration, as well as angiogenesis through the PI3K/AKT and Ras/MAPK signaling pathways [29]. In this study, GAB1 was found to be involved in several signaling pathways including the ras signaling pathway, phospholipase D signaling pathway, and ErbB signaling pathway. The expression level of GAB1 showed a positive correlation with PI(6-keto-PGF1alpha/16:0) and Cytidine-5’-diphosphocholine levels, indicating its potential importance in HCC. Additionally, previous studies have demonstrated that over expression of GAB1 was associated with poor prognosis in patients with anterior skull base meningioma, HCC and breast cancer, suggesting its potential utility as a biomarker for HCC diagnosis and prognosis [30,31,32]. Our study revealed that ME2 primarily participates in the conversion of malic acid to pyruvate and plays a crucial role in carbon metabolism. Moreover, an interaction was observed between ME2 and citric acid as well as malic acid, which are two vital components of the tricarboxylic acid (TCA) cycle. Numerous studies have also substantiated that ME2 serves as a significant source of NADPH, thereby regulating cellular reduction-oxidation (REDOX) balance and indirectly supporting fatty acids and nucleotide synthesis. Over expression of ME2 is closely associated with metabolic reprogramming in cancer, including glycolysis, amino acid metabolism, and lipid synthesis. Simultaneously, it exerts a pivotal function in the initiation, invasion, and metastasis of various malignant tumors [33]. Furthermore, specific investigations on ME2 in HCC have demonstrated its substantial involvement in glycolysis and lipid synthesis among HCC patients [34, 35]. Pyruvate kinase (PK) is a pivotal enzyme in the glycolysis pathway, catalyzing the final step by transferring a phosphate group from phosphoenolpyruvic acid (PEP) to ADP, resulting in the production of pyruvate and ATP [36]. The PKM gene encodes two isoforms, PKM1 and PKM2, with PKM2 being overexpressed in HCC and other cancers. It plays a crucial role in promoting tumor cell proliferation and metastasis, directly influencing the Warburg effect in tumor metabolism [37]. This study specifically focused on liver tissues from HCC patients and revealed high expression levels of PKM (without specifying subtypes). Furthermore, it demonstrated that PKM was involved in several key metabolic pathways including pyruvate metabolism, carbon metabolism, biosynthesis of amino acids, and glycolysis/gluconeogenesis. Metabolically linked to citric acid and malic acid - important components of the TCA cycle - through its product pyruvate. Considering PKM2 is predominantly found in tissues with high glucose metabolism such as liver tissue, and previous studies have suggested its potential use as a biomarker for malignant tumors like liver cancer [38]. Our findings further support the significance of members within the PKM family as potential biomarkers for HCC.

Recent studies have demonstrated the dynamic nature of metabolic reprogramming during carcinogenesis and cancer development. Cancer cells necessitate extensive cellular metabolic alterations to survive and proliferate under diverse environmental conditions. Consequently, metabolomics has emerged as a prominent research avenue in tumor cell biology for investigating metabolic reprogramming in tumors [39]. As a prototypical metabolism-related tumor, HCC exhibits distinct characteristics compared to normal liver tissues concerning glucose metabolism, lipid metabolism, and protein metabolism. Under normal physiological conditions, cells efficiently utilize glucose as an energy source through orderly metabolic processes such as glycolysis, the TCA cycle, and oxidative phosphorylation to generate ATP supporting various cellular activities. However, HCC displays significant changes in glucose metabolism. HCC cells exhibit enhanced glucose uptake capacity leading to heightened glycolytic activity and increased lactate production even in oxygen-rich environments. This results in an elevated ratio of lactate to pyruvate—an abnormal phenomenon known as the “Warburg effect“ [37]. Furthermore, the abundant metabolic intermediates generated during glycolysis not only fulfill immediate energy demands but also provide essential biosynthetic precursors facilitating rapid proliferation and division of HCC cells. The comprehensive analysis combining proteomics and metabolomics conducted in this study further elucidated the intricate nature of energy metabolism in HCC. Among all detected metabolic pathways, pyruvate metabolism pathway holds significant importance for energy production in HCC with its enrichment ranking among the top 10 pathways identified. It is noteworthy that citric acid serves as a pivotal intersection point connecting multiple energy metabolism pathways [40, 41].

In this study, PI (6 keto-PGF1alpha/16:0), 13,16, 19-docosatrienoic acid, PA (18:2(9Z,12Z)/20:1(11Z)) and PG (20:3(6,8,11)-OH(5)/18:2(9Z,12Z)) exhibited significant up-regulation in HCC. They demonstrated a strong diagnostic ability for HCC (AUC > 0.90) and were found to be correlated with several key node proteins such as B4GALT5. Although no studies have been conducted on these four molecules specifically, they are likely to play a pivotal role in the rearrangement of lipid metabolism in HCC as representatives of different subclasses within the lipid metabolism network. Furthermore, multiple metabolomics studies on HCC have indicated an up-regulation of phosphatidylinositol, phosphatidic acid, glycerophospholipid, polyunsaturated fatty acids and their derivatives [12, 42, 43, 44]. These abnormal alterations in lipid metabolism may be closely associated with liver cancer cell proliferation and tumor angiogenesis. Reprogramming of amino acid metabolism is a crucial characteristic of tumor metabolic abnormalities. Research has revealed that amino acids not only serve as nitrogen donors to provide raw materials for the synthesis of biological macromolecules, such as nucleotides required for tumor cell proliferation, invasion, and immune evasion but also act as vital metabolic substances for immune cells to activate and exert an anti-tumor role within the tumor microenvironment. Abnormal alterations in amino acid metabolism are closely associated with the occurrence and progression of tumors and tumor immunity. Furthermore, certain key proteins or enzymes involved in these metabolic pathways can be utilized as biomarkers for tumor diagnosis and treatment [45, 46]. Among the top 10 metabolic pathways identified in KEGG enrichment analysis conducted in this study, 4 were associated with amino acid metabolism. These included tryptophan metabolism, alanine, aspartate and glutamate metabolism, arginine biosynthesis, and arginine and proline metabolism. Previous studies have also demonstrated the crucial role of these four metabolic pathways in HCC, specifically in cell proliferation, energy supply, and immune evasion. Notably, impaired tryptophan metabolism can hinder NAD + synthesis in the liver and promote liver tumorigenesis through DNA damage [47]. The metabolism of alanine, aspartate, and glutamate is not only associated with changes in early HCC patient characteristics but also linked to cell proliferation and intracellular signal transduction in HCC [45, 48]. Dirk et al. [49] demonstrated that arginine biosynthesis was enhanced in HCC patients, accompanied by RNA binding motif protein 39 (RBM39) regulation on metabolic-related genes, promoting the reprogramming of metabolic pathways in HCC. In this study, N2-acetylornithine was considered as a potential biomarker for HCC due to its critical role in arginine biosynthesis. The metabolism of arginine and proline is not only closely associated with ATP production, protein synthesis, and nucleotide synthesis in tumor cells, but also plays a role in maintaining intracellular REDOX homeostasis [50, 51]. The degradation of arginine during the urea cycle produces several metabolites including urea, ornithine, and polyamines. Polyamines are a class of small organic polycations that are crucial for DNA replication, protein translation, and cell proliferation [52]. Among them, spermidine-a key member of the polyamine family—exerts essential biochemical characteristics in controlling cell growth, regulating the cell cycle, providing antioxidant protection mechanisms, facilitating cell signal transduction and gene expression regulation [53]. Therefore, spermidine is considered as a potential biomarker for HCC in this study.

This study conducted a comprehensive protein-metabolite analysis to explore the molecular characteristics of HCC. Our findings indicated the presence of abnormal changes in 10 pathways, including tryptophan metabolism, and suggested the potential of 12 biomarkers as candidates for HCC diagnosis or prognosis. These results suggested novel avenues for understanding the pathogenesis of HCC and could potentially contribute to the discovery of HCC-specific biomarkers. However, it is important to acknowledge certain limitations in our study. Firstly, this study was conducted at a single center with a relatively small sample size. The limited sample size may somewhat impact the generalizability of our findings. Secondly, we did not integrate the omics results with the clinical characteristics of HCC patients in this study. Further research should aim to incorporate clinicopathological features and patient prognosis data for a more comprehensive association analysis and evaluation of the clinical value of these biomarkers. Thirdly, the results of this study have not been validated in other cohorts, and extensive prospective cohort studies are required to further validate our findings.

Conclusions

The present study presented a comprehensive analysis of the intricate alterations in HCC at the proteome and metabolome levels through an integrated multi-omics investigation. We have identified 10 aberrant metabolic pathways and discovered 12 potential diagnostic markers, thereby providing novel insights into the pathogenesis of HCC and offering promising candidates for HCC biomarkers.

Table 3 Differentially expressed metabolites involved in the top 10 enrichment pathways

Functional analysis of DEMs

In order to determine the specific functions and roles of DEMs in biological processes, this study performed pathway annotation and enrichment analysis of 500 DEMs based on KEGG database. KEGG annotation results showed that 41% of the annotated pathways belonged to the metabolic category, in which lipid metabolism, amino acid metabolism, cofactor and vitamin metabolism, and carbohydrate metabolism had the largest number of metabolites (Fig. 6A). KEGG enrichment analysis showed that the DEMs were mainly enriched in glycerophospholipid metabolism, choline metabolism in cancer, phospholipase D signaling pathway, and aldosterone synthesis and secretion pathways (Fig. 5B). Among the top 20 pathways exhibiting the lowest adjusted p-value, the overall expression of DEMs involved in 8 specific pathways, including glycerophospholipid metabolism, choline metabolism in cancer, phospholipase D signaling pathway, and arginine and proline metabolism, tended to be up-regulated. Conversely, the overall expression of DEMs involved in another 12 pathways, such as aldosterone synthesis and secretion, arachidonic acid metabolism, arginine biosynthesis, and oxidative phosphorylation, exhibited a tendency for down-regulation (Fig. 6C).

Fig. 6

Functional analysis of differentially expressed metabolites. (A) KEGG annotation analysis was performed to identify differentially expressed metabolites. (B) KEGG enrichment analysis was carried out to determine the significance of these differentially expressed metabolites in metabolic pathways. Abbreviations: M, metabolism; HD, human diseases; EIP, environmental information processing; OS, organismal systems. (C) Differential abundance score map of KEGG pathways. The horizontal axis represents the DA Score and the vertical axis indicates the name of each metabolic pathway in KEGG database. The DA Score reflects overall changes in all annotated metabolites within a given pathway, with a score of 1 indicating up-regulation and − 1 indicating down-regulation for all annotated differential metabolites within that pathway. The length of line segments corresponds to absolute values of DA Scores, while dot sizes indicate numbers of differential metabolites annotated in each pathway

Integrated analysis of proteomics and metabolomics

Spearman correlation analysis

In this study, both proteomics and metabolomics analyses were conducted on HCT and DNT samples from six HCC patients, and the resulting datasets were merged for a comprehensive two-omics analysis. The two-way orthogonal partial least squares (O2PLS) method was employed to assess the intrinsic correlation between the two datasets. The combined score plot shows great similarity between the two data matrices (Fig. 7A). To gain a more intuitive understanding of the co-regulatory relationship between DEPs and DEMs, spearman correlation coefficient calculations were performed between each DEP and DEM, resulting in a correlation coefficient matrix. Heat maps were then generated to visualize the correlation strengths, focusing on the top 20 DEPs from the proteomics dataset and top 20 DEMs from the metabolomics dataset, respectively, based on their smallest adjusted p-values. To gain a more intuitive understanding of the co-regulatory relationship between DEPs and DEMs, a correlation coefficient matrix was derived by computing the spearman correlation coefficients between each pair of DEPs and DEMs. This matrix was then visualized using heatmaps, which displayed the correlation strengths for the top 20 DEPs and top 20 DEMs with the smallest adjusted p-values from their respective omics datasets (Fig. 7B). Furthermore, within this matrix, the top 20 protein-metabolite pairs with the highest absolute correlation coefficients were selected, and their correlations were plotted using a chord diagram to illustrate the strong interactions (Fig. 7C). The results showed that B4GALT5 was positively correlated with four metabolites, including PG (18:2(9Z,12Z)/18:1(9Z)), PG (20:3(6,8,11)-OH(5)/18:2(9Z,12Z)), PI (6 keto-PGF1alpha/16:0), and PA (18:2(9Z,12Z)/20:1(11Z)), while LRP4 was negatively correlated with 5,8,11-Eicosatrienoic acid, Tetracosapentenoic acid, Citric Acid, and 13,16,19-Docosatrienoic acid.

Combined analysis of biological functions

Based on the KEGG database, this study compared the pathways involved by DEPs and DEMs and found that they jointly participate in 178 pathways (Fig. 7D). The KEGG annotation results indicated that the most significant pathways, considering the number of proteins and metabolites involved, include Bile secretion, Biosynthesis of amino acids, Glycerophospholipid metabolism, and Tryptophan metabolism (Fig. 7E). The KEGG enrichment analysis revealed that the top 10 pathways exhibiting significant enrichment including tryptophan metabolism, butanoate metabolism, pyruvate metabolism, and bile secretion. Notably, both tryptophan metabolism and bile secretion pathways demonstrated a high number of DEPs and DEMs as well as substantial degrees of enrichment (Fig. 7F, Additional file 1 and Table 3). Subsequently, this study investigated the interaction between DEPs and DEMs in the top 10 enriched pathways. The results revealed that L-Glutamic acid, L-Ornithine, Malic acid, and Citric acid exhibited interactions with GLUD1 and other proteins. Notably, Citric acid demonstrated an up-regulation trend in HCC, while the other three metabolites displayed a down-regulation trend in HCC. Among these differential metabolite-interacting proteins, GLS2 showed the highest level of interaction and was found to be down-regulated in HCC (Fig. 8A and B).

Fig. 7

Consistency analysis of differentially expressed proteins and metabolites. (A) Combined sample score plot of the O2PLS model for the protein-metabolism-set. The squares and dots represented samples from different groups. The colors of the dots indicated different omics samples, with blue representing protein samples (denoted by ‘t’ after the sample name) and green representing metabolome samples (denoted by ‘u’ after the sample name). In this figure, the horizontal and vertical coordinates represented the combined score of the metabolome and proteome. (B) Correlation analysis on the top 20 differentially expressed proteins and metabolites between HCT and DNT. r represents the correlation value. (C) Chord diagram of the top 20 protein-metabolite pairs ranked by their correlation scores. Red represents positive correlation and blue represents negative correlation. (D) Venn diagram of differentially expressed proteins and metabolites involved in the KEGG pathway. (E) KEGG annotation analysis of differentially expressed proteins and metabolites. (F) KEGG enrichment analysis of differentially expressed proteins and metabolites

Fig. 8

Interaction analysis of the differentially expressed proteins and metabolites involved in the top 10 enrichment pathways. (A) Distribution of the differentially expressed proteins and metabolites involved in the top 10 enrichment pathways. (B) The interaction between the differentially expressed proteins and metabolites involved in the top 10 enrichment pathways. The diamonds represent metabolites, while the circles represent proteins. The size and color of the figures indicate the number of interacting molecules with them, with larger and darker figures representing more interacting molecules

Screening of potential biomarkers

To further identify potential biomarkers highly associated with the occurrence of HCC, the following screening criteria were formulated in this study: (1) The DEPs (DEMs) between HCT and DNT exhibited statistically significant differences (adjusted p-value < 0.05) and demonstrated a trend of high expression in HCT. (2) The absolute correlation coefficient between the DEPs and DEMs was greater than 0.85. (3) These biomarkers participated in the top 10 enriched pathways, and their interaction degree was greater than or equal to 2. (4) The area under curve (AUC) of these molecules for distinguishing HCC was greater than or equal to 0.90. Based on the above screening conditions and literature review results, a total of twelve potential biomarkers for HCC were identified in this study, including five proteins (PTP4A3, B4GALT5, GAB1, ME2, and PKM) and seven metabolites (PI(6 keto-PGF1alpha/16:0), 13,16,19-docosatrienoic acid, PA(18:2(9Z,12Z)/20:1(11Z)), Citric Acid, PG(20:3(6, 8, 11)-OH(5)/18: 2(9Z, 12Z)), Spermidine, and N2-Acetylornithine) (Table 4 and Additional file 2).

Table 4 The diagnostic value of potential biomarkers