GSE63067, GSE89632, and GSE102079 datasets were downloaded from Gene Expression Omnibus (GEO), an open-access database that provides gene expression profiles. GSE63067 (Frades et al., 2015) and GSE102079 (N. Chiyonobu, S. Shimada, Y. Akiyama, K. Mogushi, M. Itoh, K. Akahoshi, S. Matsumura, K. Ogawa, H. Ono, Y. Mitsunori, D. Ban, A. Kudo, S. Arii, T. Suganami, S. Yamaoka, Y. Ogawa, M. Tanabe and S. Tanaka, 2018) are both based on the GPL570 [(HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array]. GSE89632 (B. M. Arendt et al., 2015) is based on [(GPL14951) Illumina HumanHT-12 WG-DASL V4.0 R2 expression bead chip]. The title of the GSE63067 data set is “Expression data from human non-alcoholic fatty liver disease stages”. The data contained the gene expression profiles of 11 NAFLD patients and seven non-NAFLD controls. The title of the GSE102079 data set is “FABP4 overexpressed in intratumoral hepatic stellate cells within hepatocellular carcinoma with metabolic risk factors”. Between 2006 and 2011, 152 patients who underwent curative hepatic resection for HCC at Tokyo Medical and Dental University Hospital participated in an integrated gene expression microarray study. In the control group, 14 adjacent liver tissues were obtained from patients with metastases of colorectal cancer without chemotherapy. The validation data set was from GSE89632 and The Cancer Genome Atlas (TCGA) data set. The title of the GSE89632 data set is “Genome-wide analysis of hepatic gene expression in patients with non-alcoholic fatty liver disease and healthy donors with hepatic fatty acid composition and other nutritional factors”. A cross-sectional study included 20 patients with simple steatosis (SS), 19 non-alcoholic steatohepatitis (NASH), and 24 healthy liver donors. The TCGA database of liver hepatocellular carcinoma (LIHC) contains RNA-Seq data for 374 HCC patients and 50 normal tissues (https://portal.gdc.cancer.gov/) for gene expression and immune system infiltrates.

A comparison of DEGs between NAFLD and normal controls, HCC patients with MRFs, and corresponding controls was performed using the limma R package“complexheatmap” and “ggplot2” to generate heat maps and volcano maps, respectively, which is an efficient analysis method in bioinformatics (M. E. Ritchie et al., 2015). In NAFLD datasets, the selected criteria were p-value <0.05 and |log2FC|>1. In HCC datasets, the selected criteria were p-value <0.05 and |log2FC|>1. Additionally, the overlapping DEGs between NAFLD and HCC with MRFs were determined by Venn diagrams using the Venn platform (http://bioinformatics.psb.ugent.be/webtools/Venn/). A subsequent analysis was performed on these overlapping DEGs.

GO function enrichment analyses were conducted on the above overlapping DEGs. It consisted of biological process (BP), cellular component (CC), and molecular function (MF) (2006). The analysis of KEGG signaling pathway enrichment using a package called “clusterProfiler” (M. Kanehisa et al., 2016). GO terms and KEGG pathways enriched with adjusted p-value of 0.05 were selected for analysis.

In order to further investigate the interactions between the above-mentioned common genes, a search tool called the Search Tool for Retrieval of Interacting Genes (STRING) has been developed for PPI network construction (D. Szklarczyk et al., 2015). Interaction scores of at least 0.4 were considered significant. Subsequently, PPI network visualisation was conducted with Cytoscape software. Then, the Maximal Clique Centrality (MCC), Density of Maximum Neighborhood Component (DMNC), Maximum Neighborhood Component (MNC), Degree, and Edge Percolated Component (EPC), algorithms in the cytoHubba plug-in (http://hub.iis.sinica.edu.tw/cytohubba/) was applied to identify PPI hub genes with high connectivity.

Based on the TCGA database, the six hub genes expression in HCC normal tissues and tumor tissues was investigated. In LIHC, 374 HCC specimens with normal adjacent tissues and HCC tissue (50 each) were compared with neighbouring normal tissues. GEPIA was used to investigate the prognostic significance of hub genes (http://gepia.cancer-pku.cn/index.html) (Z. Tang et al., 2017). Survival analyses were considered significant when log-rank p < 0.05 was used.

Based on the training dataset, hub genes associated with prognosis were identified and assessed against other datasets for their predictive performance. Half of TCGA is set as the training set. The other half of TCGA is set as a validation set. The entire TCGA cohort is a verification set. Using univariate Cox proportional hazard regression analysis, it was evaluated whether hub genes are associated with overall survival (OS) in the training set. In the “glmnet” package, the Latent Selection Operator penalised Cox proportional hazard regression using Cox proportional hazards models. A prediction formula for gene characteristics was devised. The formula for the model is as follows: risk score = gene1×β1 (gene one expression level) + gene 2×β2 (gene two expression level) +…gene n×βn (gene n expression level). In this formula, genes are combined with gene expression values and regression coefficients from multiple Cox proportional hazards regression models (George et al., 2014; Zhang J. et al., 2021). Using the Kaplan–Meier (K–M) survival curves, survival comparisons were performed between low- and high-risk groups via the R package “survival”. Furthermore, a time-dependent receiver operating characteristic (ROC) analysis (including 1-, 3-, and 5-year survival) was conducted to evaluate hub gene sensitivity and specificity using the R package “survival ROC” (P. J. Heagerty et al., 2000). It is critical to consider the area of the AUC curve when trying to predict clinical outcomes. Prognosis is better when AUC >0.5; the closer AUC is to 1, the better.

Tumor contains a large number of immune cells, and the prognosis of high-grade HCC patients with high subtype of dominant immunity is obviously better (Y. Kurebayashi et al., 2018). To examine whether the expression of hub genes is correlated with the presence of immune cells in HCC, we examined the correlation between hub gene mRNA expression and tumor-infiltrating immune cells. The web tool TIMER was used (https://cistrome.shinyapps.io/timer/) (T. Li et al., 2017). Six tumor-infiltrating cell subsets were analysed, such as B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils, and dendritic cells.

ROC curves for hub genes were constructed using the pROC package (X. Robin et al., 2011). To measure the effectiveness of the model, we calculated the area under the curves (AUC). These results showed the usefulness of genes for diagnostic purposes. p-values <0.05 were considered statistically significant.

The series GSE63067 dataset about NAFLD and the series GSE102079 dataset about HCC from the NCBI GEO database was downloaded. Based on a p-value of 0.05 and |log2FC | of >1.0, 125 DEGs were identified in GSE63067, and 726 DEGs were identified in GSE102079 using the “limma” package in R software. Volcano plots and heatmaps were used to visualise the DEGs of the two data sets shown in Figures 1A,B and Figures 2A,B, respectively. Using the Venn Diagram online tool, 26 common genes of two diseases were identified and are shown in Figure 1C.

Functional enrichment and KEGG pathway analyses of 26 common NAFLD and HCC genes were performed at a threshold of p-value <0.05. The results showed that DEGs were enriched in biological processes, including cellular response to environmental stimuli, cellular response to abiotic stimuli, cellular response to ionising radiation, unsaturated fatty acid biosynthetic process, and response to zinc ions (Figure 3A). Regarding molecular function, DEGs were principally associated with receptor ligand activity, cytokine activity, monocarboxylic acid binding, and fatty acid binding Figure 3B. The KEGG pathways of DEGs were enriched in the PPAR signalling pathway (Figure 3C).

The PPI network was first performed based on the STRING database to investigate how DEGs interact with one another. Afterwards, the results were imported into Cytoscape to be analysed (Figure 4A). A Cytoscape plug-in, Cytohubba, was used to analyse the PPI network and identify hub genes. We got top 10 genes from protein-protein network ranked by five different algorithms of cytohubba including MCC, DMNC, MNC, Degree, and EPC (Table1). In this study, the genes with the top six values were considered as hub genes. Based on the five algorithms, the top six genes were determined to be hub genes: FABP5, SCD, CCL20, AGPAT9, PLIN1, IL1RN (Figure 4B).

To evaluate the diagnostic value of the top six hub genes obtained from the above analysis, ROC curves were constructed and their corresponding area under the curve (AUC) was calculated. Figure 5A shows the result of NAFLD. The AUC for FABP5, SCD, CCL20, AGPAT9, PLIN1, and IL1RN in NAFLD patients and normal controls were 0.828, 0.818, 0.883, 0.857, 0.961, 0.818 at NAFLD GSE63067 dataset. Figure 5B shows the ROC curves in HCC patients and normal controls. The AUC for FABP5, SCD, CCL20, AGPAT9, PLIN1, and IL1RN in HCC and the normal controls were 0.717, 0.785, 0.776, 0.821, 0.901, 0.839 at HCC with MRFs GSE102079 dataset. For validation, The AUC for FABP5, SCD, CCL20, AGPAT9, and IL1RN were 0.561, 0.746, 0.858, 0.981, 0.907 in NAFLD based on GSE89632 (Figure 5C). The AUC for FABP5, SCD, CCL20, AGPAT9, PLIN1, and IL1RN in HCC and the normal controls were 0.900, 0.614, 0.807, 0.589, 0.725, 0.769 at HCC TGCA dataset Figure 5D).

Consistent with the results in GEO datasets, the mRNA expression of AGPAT9, PLIN1, and IL1RN was significantly downregulated. At the same time, that of FABP5, SCD, and CCL20 were upregulated considerably in TCGA HCC compared with non-tumor tissues (Figures 6A–L). According to the GEPIA web tool, FABP5 and PLIN1 mRNA expression are significantly linked to overall survival (OS) (Figures 6M–Q).

The correlation of expression levels of hub genes was captured using GEPIA. An analysis of correlation was performed on any two of FABP5, SCD, CCL20, AGPAT9(GPAT3), PLIN1, and IL1RN and six hub genes. The above data indicate that upregulation of one of them will decrease the high expression of other genes. IL1RN and GPAT3 (Figure 7A), CCL20 and SCD (Figure 7B), CCL20 and FABP5 (Figure 7C), SCD and PLIN1 (Figure 7D), IL1RN and CCL20 (Figure 7E), and IL1RN and PLIN1 (Figure 7F) are all positively related to each other. This may indicate that there is a common transcription factor as well as epigenetic modifications controlling them all.

To avoid overfitting by LASSO regression, lambda. Min was selected, resulting in a more accurate prediction rate. We used the multivariate Cox proportional hazards regression analysis. Five prognostic genes were developed, including ABP5, CCL20, GPAT3, PLIN1, and IL1RNIn order to calculate the risk score for each patient, the following formula was used: risk score = (0.170124209 ×FABP5) + (0.073621309 ×CCL20) + (0.011005683 × GPAT3) + (−0.056212587 ×PLIN1)+ (−0.100588077 ×IL1RN). The LASSO coefficient for SCD is equal to 0. Hub gene risk scores were used to determine whether HCC patients were low-risk or high-risk (Figure 8A). A significantly worse OS was observed in high-risk patients compared to low-risk patients (Figure 8B, training set p = 0.008, validation set p = 0.0026, entire TCGA set p = 0.001). The reliability of hub genes was subsequently assessed using time-dependent ROC curves (Figure 8C). As a result, the area under the curve (AUC) was 0.738, 0.612, and 0.695 for 1-year, 3-year, and 5-year survival, respectively for the training set. The AUC was 0.611, 0.633, and 0.664 for 1-year, 3-year, and 5-year survival, respectively, for the validation set. These curves were also applied in the entire TCGA set. The AUC was 0.696, 0.634, and 0.673 for 1-year, 3-year, and 5-year survival, respectively.

There are three kinds of cells in the tumor microenvironment: tumor cells, stromal cells, and immune cells that infiltrate the tumor. The TIMER web tool showed that the expression of all six hub genes was associated with infiltrating immune subsets, and the expression of NAFLD and HCC showed the most significant correlation with them. For B-cells, CD4⁺ T-cells, CD8⁺ T-cells, neutrophils, macrophages, and dendritic cells, the expression of FABP5 showed the most significant correlation with them (Figures 9A–F).