With profound modern lifestyle changes, non-alcoholic fatty liver disease (NAFLD) has become an epidemic, imposing a heavy burden on global public health [1]. The non-alcoholic fatty liver disease encompasses a spectrum of liver disorders from simple hepatic steatosis to non-alcoholic steatohepatitis (NASH), and NASH patients are at high risk of cirrhosis and hepatocellular carcinoma (HCC) [2],[3]. Hepatic fibrosis, though not a diagnostic criterion for NAFLD, is frequently observed and linked to poor clinical outcomes [4].

Identifying disease-specific gene expression patterns elucidates NAFLD molecular mechanisms. It enables development of reliable diagnostic and therapeutic response biomarkers. Tools including WGCNA analyze high-dimensional datasets. The weighted gene co-expression network analysis screens therapeutic targets and diagnostic biomarkers [5],[6],[7]. Mendelian randomization (MR) is a robust method for inferring causal relationships between exposure factors and disease outcomes, overcoming confounding factors and reverse causality bias via single nucleotide polymorphisms (SNPs) as instrumental variables (IVs) [8],[9].

In the present research, we conducted an in-depth analysis of liver tissue transcriptomic profiles derived from NAFLD patients and healthy controls (HC). We first screened out DEGs between the two cohorts and then applied weighted WGCNA to mine gene modules that are closely associated with NAFLD pathogenesis. On this basis, we pinpointed hub genes with potential diagnostic value for NAFLD and further explored the causal link between FOS (Fos proto-oncogene) expression and NAFLD susceptibility through MR analysis. Our work is intended to yield novel insights into the clinical diagnosis of NAFLD and the in-depth investigation of its underlying molecular mechanisms.

Data source

We retrieved the gene expression profile dataset GSE89632 from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), which encompasses detailed clinical characteristics as well as liver tissue gene expression data obtained via transcriptome microarray assays from 19 patients diagnosed with NAFLD and 24 healthy control subjects. All raw gene expression data were obtained directly from the GEO database.

Differentially expressed genes identification

Initially, data from the GSE89632 dataset were read using R software (version 4.3.0) and then normalized. Subsequently, DEGs were screened with the “limma” package, applying a differential threshold of |log2FC| > 1 and FDR < 0.05. After statistical analysis to determine significance, volcano plots and heatmaps of DEGs were generated using the “pheatmap” and “ggplot2” packages in R (Figure 1).

Weighted gene co-expression network analysis

The R package “WGCNA” was used to construct a NAFLD-related gene co-expression network. A soft-thresholding power of 6 was selected based on the scale-free topology criterion (R² = 0.88), which ensured the constructed gene co-expression network approximated a scale-free network structure and minimized spurious correlations between genes. The correlation between different gene modules and NAFLD pathogenesis was evaluated, and the most disease-relevant module was selected to extract candidate hub genes (Figure 2). By quantifying pairwise gene expression correlations and their associations with phenotypic traits, WGCNA effectively identifies co-expressed gene modules with coordinated biological functions, making it a robust tool for prioritizing potential disease biomarkers.

Identification of candidate pivotal genes and GO/KEGG analysis

Genes overlapping between WGCNA-derived module genes and DEGs were defined as potential NAFLD-associated hub genes. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the overlapping genes were performed using the R package “clusterProfiler” to explore NAFLD progression mechanisms, with the KEGG database used for systematic gene function analysis (Figure 3A, Figure 3B, Figure 3C).

Discovery of hub genes in the protein-protein interaction network

The STRING database (version 2) was used to construct protein–protein interaction (PPI) networks, which were visualized via Cytoscape software. The Degree algorithm in Cytoscape was applied to rank genes by connectivity in the PPI network, and the top-ranked genes were identified as NAFLD hub genes (Figure 4).

Nomogram model construction

A nomogram model for NAFLD risk prediction was constructed using the R package “rms”. Harrell’s concordance index (C-index) was calculated to evaluate predictive power. The R package “pROC” was used to plot ROC curves and calculate the area under the ROC curve (AUC), with 0.9 ≤ AUC < 1 defined as excellent diagnostic accuracy (Figure 5).

Immune cell infiltration analysis

We adopted the CIBERSORT algorithm to quantify the relative proportions of 22 distinct immune cell subpopulations, thereby characterizing the landscape of immune cell infiltration in NAFLD tissue samples. Furthermore, we performed a correlation analysis between FOS gene expression levels and the infiltration abundance of each immune cell subtype. Patients were divided into FOS high- and low-expression groups by median split, aiming to clarify the potential immunomodulatory function of FOS in the pathological progression of NAFLD (Figure 6).

Mendelian randomization analysis

Two-sample MR analyses were conducted with the R package “TwoSampleMR” to probe into the causal relationship between FOS gene expression levels and the risk of developing NAFLD, in which SNPs were employed as IVs. Genome-wide association study (GWAS) data of FOS were retrieved from public repositories. The inverse variance weighting (IVW) method was the primary analysis, and MR-Egger regression was used for sensitivity analysis to detect horizontal pleiotropy. Leave-one-out sensitivity analysis was also conducted to verify result robustness (Figure 7). All data used were obtained from public databases.

Screening for DEGs

We screened DEGs between NAFLD and control groups in GSE89632. We identified 387 DEGs: 162 upregulated and 225 downregulated in NAFLD samples (Figure 1).

Construction of WGCNA network and identification of NAFLD-related module

We applied WGCNA to the GSE89632 dataset. We identified 12 independent gene co-expression modules with distinct expression patterns [10]. We calculated correlation coefficients to evaluate module–NAFLD associations. The cyan module showed the strongest positive correlation with NAFLD. We defined the cyan module as the NAFLD-related core module. It regulates NAFLD initiation and progression (Figure 2).

GO and KEGG enrichment analyses

A total of 115 overlapping genes were screened out from the intersection of the cyan module genes identified via WGCNA and the DEGs, and these genes were defined as potential hub genes associated with NAFLD (Figure 3A). The GO enrichment analysis showed these genes were mainly involved in biological processes such as fat cell differentiation, amoeboid-type cell migration, and response to lipopolysaccharide/bacterial origin molecules (Figure 3B). The Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis indicated enrichment in the TNF signaling pathway, MAPK signaling pathway, AGE-RAGE signaling pathway in diabetic complications, osteoclast differentiation, and IL-17 signaling pathway (Figure 3C) [11].

PPI network analysis for hub genes

The STRING online tool constructed a PPI network of the overlapping hub genes, and Cytoscape identified the top 10 upregulated genes with the highest connectivity scores: FOS, JUN, NR4A1, JUNB, EGR1, MYC, IL1B, CCL2, CXCL8, and PTGS2 (Figure 4). Color intensity in the plot represents connectivity score, with deeper colors indicating higher scores.

Development of a nomogram model for NAFLD risk stratification

A nomogram model incorporating the identified hub genes was established to predict the risk of NAFLD, and this model demonstrated robust predictive performance in clinical sample analysis (Figure 5A). We further conducted ROC curve analysis to evaluate the individual diagnostic value of five key hub genes (FOS, JUN, IL1B, MYC, EGR1), which all exhibited high diagnostic efficacy: the area under the ROC curve (AUC) values were 1.000 for FOS, 1.000 for JUN, 0.939 for IL1B, 1.000 for MYC, and 0.991 for EGR1, respectively (Figure 5B). The AUC value of the constructed nomogram also verified its strong capability to discriminate NAFLD patients from healthy control subjects accurately.

Analysis of immune cell infiltration characteristics in NAFLD tissues

Immune cell infiltration profiling uncovered distinct immune landscape patterns that were closely associated with FOS expression levels in NAFLD tissue samples (Figure 6). Samples were divided into FOS high-expression and low-expression groups using the median expression level of FOS as the cutoff value. In comparison with the FOS low-expression group, the FOS high-expression group presented markedly elevated infiltration abundances of γδ T cells, M1 macrophages, M2 macrophages, resting dendritic cells and resting mast cells. Conversely, the infiltration levels of naive B cells, plasma cells, monocytes, activated dendritic cells, activated mast cells and neutrophils were significantly lower in the FOS high-expression group. These findings confirm a direct and close correlation between FOS gene expression and the immune cell infiltration profile in the NAFLD microenvironment.

No causal association between FOS and NAFLD risk

The F-statistic of instrumental variables was greater than 10, indicating no weak instrument bias. Inverse variance weighting analysis showed no significant causal association between FOS levels and NAFLD risk (OR = 0.997, 95% CI = 0.947–1.049, p = 0.898) (Figure 7A and Figure 7B). The Mendelian randomization-Egger regression also showed no statistically significant causal relationship (OR = 1.021, 95% CI = 0.891–1.170, p = 0.774). The funnel plot showed a roughly symmetrical distribution of SNP effect sizes (Figure 7C), and the MR-Egger regression intercept was not statistically significant (p = 0.719), indicating no horizontal pleiotropy. Leave-one-out sensitivity analysis yielded consistent outcomes across all iterations (Figure 7D), which further validated the robustness and reliability of the Mendelian randomization analysis results.

The non-alcoholic fatty liver disease is a chronic liver disorder. It presents as excessive lipid accumulation in hepatocytes without heavy alcohol intake. Its molecular mechanisms remain unclear.

The non-alcoholic fatty liver disease progresses sequentially: fibrosis, cirrhosis, and HCC in severe cases. Its pathogenesis involves genetic susceptibility, metabolic dysregulation, and environmental factors [1],[2]. Lifestyle modification is first-line therapy. Targeted molecular treatments are urgently needed [12]. Molecular biology advances drive NAFLD gene research. Research aims to map regulatory networks and identify therapeutic targets [12]. Molecular biology advances drive NAFLD gene research. Research maps regulatory networks and identifies therapeutic targets.

We applied WGCNA and identified 10 candidate hub genes (FOS, JUN, NR4A1, JUNB, EGR1, MYC, IL1B, CCL2, CXCL8, PTGS2). The hub-gene-based nomogram model predicts NAFLD susceptibility accurately. Receiver Operating Characteristic analysis confirms high diagnostic efficiency for five core hub genes (FOS, JUN, IL1B, MYC, EGR1). These findings underscore the great potential of these genes to serve as promising tools for the clinical diagnosis and clinical management of NAFLD.

Cytoscape-based analysis revealed that FOS was the top-ranked hub gene within the constructed PPI network, which implies that this gene may exert a central regulatory function in the molecular regulatory networks linked to NAFLD pathogenesis. As a classic proto-oncogene, FOS (c-fos) encodes a transcription factor belonging to the activator protein-1 (AP-1) family; this factor forms a heterodimer with Jun family proteins, thereby mediating the transcriptional regulation of downstream target genes [13],[14]. FOS can be rapidly induced by a variety of extracellular stimulatory signals and has been the subject of extensive investigation in oncology research, with its dysregulated expression being detected in a wide range of malignant tumors, including hepatic carcinoma [15],[16],[17]. In addition, FOS is involved in immune responses, inflammation, and metabolic diseases [18],[19].

A previous study demonstrated that the AP-1 complex plays a crucial role in NAFLD development: FOS promotes PPARγ expression, and different AP-1 dimers uniquely regulate PPARγ at the transcriptional level, making AP-1 a key link between obesity, hepatic lipid metabolism, and NAFLD [20]. To the best of our knowledge, this is the first study to investigate the causal relationship between FOS levels and NAFLD risk using two-sample MR analysis based on comprehensive GWAS data. Mendelian randomization results did not support a causal association between elevated FOS levels and increased NAFLD risk, with MR effectively reducing systematic biases in traditional observational studies and accurate genotyping minimizing regression dilution bias, ensuring result reliability [8],[9].

This study is subject to several limitations. Firstly, owing to the scarcity of available microarray data for NAFLD, our analysis was restricted to the single GSE89632 dataset; the integration of multiple independent datasets would undoubtedly enhance the robustness and generalizability of our findings. Secondly, hub gene identification and analysis relied solely on bioinformatics. In vitro and in vivo experiments are required to clarify their regulatory mechanisms in NAFLD. Thirdly, MR analysis used GWAS data from European populations only. Results do not generalize to other ethnic groups. Fourth, the AUC values of 1.000 observed for FOS, JUN, and MYC may reflect potential overfitting caused by the relatively small sample size (19 NAFLD vs 24 controls), which should be interpreted with caution.

We built a WGCNA-based gene co-expression network. We screened NAFLD-associated hub genes. These genes clarify NAFLD molecular regulatory pathways. Among these identified genes, the core hub genes including FOS, JUN, IL1B, MYC and EGR1 exhibited high diagnostic efficiency for NAFLD, and the nomogram model developed on this basis also demonstrated robust performance in predicting NAFLD risk. Two-sample MR found no causal link between FOS expression and NAFLD susceptibility. FOS correlated with immune cell infiltration in NAFLD tissues. FOS regulates the NAFLD immune microenvironment. These findings support targeted NAFLD interventions and personalized therapy. They improve early screening, optimize clinical management, and reduce global NAFLD burden.