VCI clinical and transcriptomic (mRNA) data were downloaded from GEO database. The filter conditions were set to ①” vascular cognitive impairment”; ② human. This study was derived from 10 “normal cognitive patients” and 10 “patients with VCI” from the microarray dataset GSE201482 in 10 cases. DEGs of VCI were screened according to the criteria of P-value < 0.05 and |log₂FC| ≥ 1.00. HCC clinical, transcriptomic (mRNA) data were downloaded from TCGA (https://portal.gdc.cancer.gov/) website, and their expression matrix was compiled and summarized by R. DEGs were screened by adjust P-value < 0.05 and |log₂FC| ≥ 1.0, and heat maps were plotted. DEGs were categorized into up-regulated, down-regulated, and non-statistically significant groups, and then imported into R for volcano plots. Bioconductor R software's bioconductor R package was used to normalize and calculate expression values for microarray data. Based on the screened DEGs, heat maps and cluster analyses were performed using the heatmap package. We transformed DEG P-value to –log10, grouped them according to log₂FC, and imported the processed data into R for volcano plotting. The HCC expression matrix data was downloaded from the Cancer Genome Atlas database. According to the pre-screened genes, the relative expression of the core genes in the HCC expression data was analyzed using “ggpubr” package and a relative expression box plot was drawn.

For multi-omics analysis, metabolically differentially expressed metabolites and DEGs were imported from HCC and VCI metabolomics into Metaboanalyst 5.0, while metabolomics genes associated with HCC and VCI were incorporated into Venny 2.1 software (http://bioinfogp.cnb.csic.es/tools/venny/index.html) for plotting Venn diagrams and obtaining metabolomic differential genes associated with HCC and VCI.

The HCC expression matrix data obtained above were subjected to a deconvolution algorithm using the CIBERSORT package, which allows estimation of the cellular composition of complex tissues based on normalized gene expression data, as well as quantification of specific cell types. The immune cell sorting perl script is used to sort the HCC-associated infiltrating immune cells. The limma R package in R was used to calculate expression values from the microarray data. HCC and paraneoplastic tissues were analyzed using the CIBERSORT package. With the CIBERSORT package, the composition of immune cells in each sample was further analyzed and histograms were plotted. Pheatmap was used to generate heat maps showing immune cell distribution. HCC immune cell infiltration co-expression map was plotted using corrplot to analyse interactions between immune cell populations. Finally, the expression of each immune cell was analyzed using the vioplot package in HCC tissues and paracancerous tissues. In relation to immunofunctional analysis, we used the “limma,” “GSVA,” “gseabase,” “pheatmap,” and “reshape2” packages to perform immunofunctional analysis on HCC metabolomics-related genes to identify targets for precision therapy.

Standardized HCC metabolomics data were merged with HCC clinical data, and univariate and multifactorial Cox prognostic survival analyses on HCC differential genes were performed using the R language packages survival, caret, glmnet, survminer, and survroc to plot survival curves for closely related key metabolomics genes. In addition, 370 cases were randomly divided into testing group (n = 185) and training group (n = 185). R software was used to perform risk survival analysis on the pre-merged data, plotting risk survival curves and receiver operator characteristic (ROC) curves for testing group, training group, and total group as previous studies (37–40). On the general clinical data of HCC, the survival package was again used to perform clinical statistical analysis and risk prognostic analysis. To investigate the risk factors associated with HCC, forest plots, and histograms were plotted for single- and multi-factor independent prognostic analysis. Further model validation on clinical subgroups was performed using the “survival” and “survminer” packages.

The scatterplot3d package was used for the principal component analysis of HCC metabolomics-related genes and HCC-related risk genes, and the clusterprofilergo.R package in R (https://www.r-project.org/) software and Perl language were used for metabolomics-related HCC GO analysis for differential genes. KEGG pathway enrichment analysis was performed using the clusterprofilerkegg.R package to analyze core pathway enrichment based on enrichment factor values and to investigate the biological functions and signaling pathways that may be associated with HCC.

TMB files were downloaded from the TCGA database and correlation tests between core genes and TMB were performed using functions, with correlation coefficients and P-values calculated. In addition, correlation analysis between metabolism-related genes and HCC tumor mutational load was performed using the survival and survminer software packages, and correlation coefficients and P-values were calculated.

The GSEA website was used to download the GO/KEGG annotation files for the whole transcriptome genes, and GO and KEGG enrichment analyses of the core genes were performed using the limma, org. Hs. eg. db, clusterprofiler, and enrichplot packages as previous researches (41–43). Analyses of cellular component (CC), molecular function (MF), biological process (BP), and KEGG pathway enrichment were conducted for the core genes. Based on the enrichment factor values, we analyzed the core pathway enrichment and examined the potential biological functions and signaling pathways of the HCC core genes.

“pRRophetic” is a package that can be used to predict phenotypes from gene expression data, predict drug sensitivity in external cell lines, and predict clinical data. In order to determine the drug sensitivity of each sample from TCGA database, we used the R package “pRRophetic” to obtain HCC metabolism-related genes after prescreening.

A total of 14 differentially expressed metabolites from HCC, including 2 up-regulated and 12 down-regulated metabolites, as well as 71 differentially expressed metabolites from VCI, including 32 up-regulated and 39 down-regulated metabolites, were identified by screening criteria of fold change (FC) >1.5 and P-value < 0.05. Data enrichment and metabolic pathway analysis were performed using Metabolanalyst 5.0. Ammonia recycling, glutathione metabolism, glutamate metabolism, and Beta A metabolism were enriched for differentially expressed metabolites of VCI (Figures 1A, B). As part of the metabolic signaling pathway enriched for aminoacyl-tRNA, glycine, serine, and threonine biosynthesis, valine, leucine, and isoleucine biosynthesis, glyoxylate and dicarboxylate metabolism, phenylalanine metabolism, phenylalanine, tyrosine, and tryptophan biosynthesis, etc (Figures 1A, B). Among the metabolic functions enriched for differentially expressed metabolites of HCC are the Warburg effect, gluconeogenesis, purine metabolism, phosphatidylethanolamine biosynthesis, phosphatidylcholine biosynthesis, and arginine and proline metabolism (Figures 1C, D). The main enriched metabolic signaling pathways include glycerophospholipid metabolism, citrate cycle (TCA cycle), pyruvate metabolism, purine metabolism, glycolysis/gluconeogenesis, and the alanine, aspartate, and glutamate metabolism (Figures 1C, D). Based on the above, it is evident that metabolism-related genes may play a role in HCC and VCI.

The TCGA website (https://portal.gdc.cancer.gov/) was used to download clinical and transcriptomic expression data related to HCC, and a total of 424 transcriptomic and 377 clinical datasets were obtained according to the predefined screening criteria. Differential expression data were screened based on adjusted P-values < 0.05 and |log₂FC| ≥ 1.0. A total of 882 differentially expressed mRNAs were screened in the dataset, including 487 up-regulated mRNAs and 395 down-regulated mRNAs; DEGs were screened for differential analysis based on P-value. The top 100 most significant differentially expressed mRNAs were screened based on P-values and plotted as heat maps (Figure 2A). Following differential analysis, P-values were –log10 transformed, grouped according to log2 FC (groups of up-regulated DEGs, down-regulated DEGs, and non-statistically significant DEGs), and imported into R for plotting volcanoes (Figure 2B). The GSE201482 dataset contained 343 differentially expressed mRNAs, including 179 up-regulated and 164 down-regulated mRNAs. The top 100 differentially expressed coding RNAs with the greatest significance were screened, and the heat map was created based on the P-value (Figure 2C). Further, –log10 (P-value) was transformed from differential analysis microarray data, –log10 (P-value) was grouped by log2 FC (groups of up-regulated DEGs, down-regulated DEGs, and non-statistically significant DEGs), and the processed data was imported into R to plot volcanoes (Figure 2D).

The differentially expressed metabolites and DEGs of HCC and VCI obtained in the previous stage were imported into Metaboanalyst 5.0 for multi-omics analysis, and 360 metabolic DEGs associated with HCC and 63 metabolic DEGs associated with VCI were obtained (Figures 3A, B). To plot the Venn diagram, the above metabolomics genes associated with HCC and VCI were imported into Venny 2.1 software (http://bioinfogp.cnb.csic.es/tools/venny/index.html). Accordingly, 8 mRNAs were selected, including 1 up-regulated and 7 down-regulated mRNAs, namely: NNMT, PHGDH, NR1I2, CYP2J2, PON1, APOC2, CCL2, and SOCS3 (Figure 3C). Figure 4 shows a heat map of those differentially metabolized DEGs related to HCC and VCI.

The 377 clinical cases were randomly divided into a test group and a training group, for clinical and statistical analysis. No significant differences were found between the two groups in terms of age, gender, and tumor stage (P-value > 0.05), and the data from the two groups were comparable, as shown in Supplementary Table 1. Based on a multifactorial regression model, we composed these six genes (NNMT, PHGDH, NR1I2, CYP2J2, PON1, APOC2, CCL2, and SOCS3) into a risk factor model called riskScore. In the COX survival prognostic model, survival time decreased with increasing riskScore in the test group and the training group (P-value < 0.05) (Figures 5A, B). Independent prognostic analyses by univariate and multifactorial regression showed that the tumor stage and riskScore were significantly associated with prognosis in HCC patients (Figures 5C, D). For prognosis prediction of HCC patients at 1, 3, and 5 years, the area under the curve (AUC) of the subject working characteristic (ROC) curve exceeded 0.59 (Figure 5E). Among the ROC curves for all HCC risk factors, the AUC of riskScore for metabolic genes related to HCC and VCI was the largest and >0.6 (Figure 5F), indicating a good sensitivity of the established survival prognostic model. Based on HCC and VCI-related metabolic genes, a nomogram prognostic model was constructed, and the test results showed survival rates of 0.929, 0.86, and 0.81 after 1 year (P-value < 0.05) (Figure 5G).

Principal component analysis of HCC metabolism-related genes and HCC-associated genes was performed using the scatterplot3d package (Figures 6A, B). GO and KEGG pathway enrichment analysis of those eight differentially metabolized DEGs related to HCC and VCI were done using Bioconductor package and clusterprofiler package in R language. The GO analysis of those eight differentially metabolized DEGs showed that their biological processes were mainly enriched in response to a drug or an exogenous drug catabolic process (Figures 7A–F). The cellular components that were mainly enriched included high-density lipoprotein particles and plasma lipoprotein particles. The molecular functions including phospholipase activator activity, lipase activator activity, and phospholipase binding were also enriched. The enriched KEGG pathways included TNF signaling pathway, Influenza A lipid and atherosclerosis pathway, linoleic acid metabolism, nicotinate and nicotinamide metabolism, glycine, serine, and threonine metabolism, and so on (Figures 7G, H). This indicates that the pathways and functions of these DEGs enrichment may be connected to the immune microenvironment and metabolism.

The obtained HCC expression matrix data were background corrected and normalized and expression values for the microarray data were calculated using the limma R package in R. The immune cell composition of each sample was further analyzed using the CIBERSORT package and histograms were plotted (Figure 8A). The immune cell distribution heat map was plotted using the pheatmap package (Figure 8B). The interaction of immune cell populations in HCC was then analyzed using the corrplot package and co-expression maps of immune cell infiltration in HCC were plotted (Figure 8C). Finally, the expression matrix data were analyzed using the Vioplot package to investigate the expression of each immune cell in HCC tissues as well as paracancerous tissues, and further plotted (Figure 8D).

Immunological functions of HCC metabolism-related genes were analyzed using limma, Gene Set Variation Analysis (GSVA), gseabase, pheatmap, and reshape2 packages. The immune functions were mainly focused on antigen-presenting cells (APC) co-inhibition, APC co-stimulation, cytokine-cytokine receptor interaction (CCR), checkpoint, cytolytic activity, human leukocyte antigens (HLA), inflammation-promoting, MHC_class_I, and parainflammation (Figure 9). This study reveals and visualizes immune infiltration and immune function associated with metabolic DEGs of HCC and VCI.

The HCC expression data and TMB files were imported into R, and the correlation between DEGs related to HCC and tumor mutation load was calculated by using the function, and the waterfall plots of high and low-risk groups were drawn according to the correlation results (Figures 10A, B). The differential analysis results suggested that the tumor mutation load in the low-risk group was significantly higher than that in the high-risk group (Figure 10C). Survival analysis of the tumor mutation load of HCC metabolism-related genes showed that the survival period of the low tumor mutation load group was significantly longer than that of the high tumor mutation load group (P-value < 0.05) (Figure 10D). In addition, combining tumor mutation load characteristics and metabolism-related genetic factors, the low-risk group with low tumor mutation load had the highest probability of survival while the high-risk group with low tumor mutation load had the lowest probability of survival (Figure 10E).

After obtaining the core genes PHGDH, NR1I2, and APOC2 from the previous differential expression analysis and survival prognosis analysis, we further investigated their relative expression in HCC. The transcriptome data of HCC were downloaded from the TCGA and the ggpubr package was used to analyze the relative expression of the identified core genes and to draw box expression maps (Figures 11A–C). The results revealed that PHGDH and NR1I2 genes exhibited low expression in HCC tumor tissues (P-value < 0.001) while APOC2 genes were highly expressed in HCC tumor tissues (P-value < 0.01).

The GO/KEGG annotation file downloaded from the GSEA website and the HCC tumor data file were read into R. Analytical operations were performed, and it was found that: the GO of gene PHGDH at HCC was enriched in CHROMATIN REMODELING, DNA PACKAGING, and PROTEIN DNA COMPLEX SUBUNIT ORGANIZATION functions (Figure 12A). The GO function of gene NR1I2 was enriched in ACTIVATION OF IMMUNE RESPONSE, ADAPTIVE IMMUNE RESPONSE BASED ON SOMATIC RECOMBINATION OF IMMUNE RECEPTORS BUILT FROM IMMUNOGLOBULIN SUPERFAMILY DOMAINS, and ALPHA BETA T CELL ACTIVATION (Figure 12B). The gene APOC2 was functionally enriched in CHROMATIN ASSEMBLY OR DISASSEMBLY, EPIDERMAL CELL DIFFERENTIATION, and INFLAMMATORY RESPONSE TO ANTIGENIC STIMULUS (Figure 12C); The main PHGDH gene-enriched KEGG pathways are OLFACTORY TRANSDUCTION, CIRCADIAN RHYTHM MAMMAL, GRAFT VS. HOST DISEASE signaling pathways (Figure 12D); The NR1I2 gene enriched KEGG pathways are mainly ANTIGEN PROCESSING AND PRESENTATION, CYTOKINE RECEPTOR INTERACTION, and CYTOSOLIC DNA SENSING signaling pathways (Figure 12E). The main KEGG pathways enriched by the APOC2 gene are OLFACTORY TRANSDUCTION, CYTOSOLIC DNA SENSING PATHWAY, and REGULATION OF AUTOPHAGY signaling pathway (Figure 12F).

To evaluate therapeutically effective drugs against DEGs related to HCC and VCI, we used the limma, ggpubr, prrophetic, and ggplot2 packages. Drugs with potential clinical efficacy against HCC metabolism-related genes include A-443654, A-770041, AP-24534, BI-2536, BMS-509744, CGP-60474, CGP-082996, CMK, cyclopamine, dasatinib, doxorubicin, etoposide, gemcitabine, GW843682X, HG-6-64-1, and JW-7-52-1 (Figures 13A–P).