Integrated Machine Learning and Multi-Omics Analysis Identifies Mitophagy-Related Core Genes and Mechanisms in Non-Alcoholic Fatty Liver Disease

Introduction

Non-alcoholic fatty liver disease (NAFLD) has emerged as the most common chronic liver disease worldwide, with its prevalence rising significantly alongside the epidemics of obesity, type 2 diabetes, and metabolic syndrome.¹ The pathological spectrum of NAFLD ranges from simple steatosis to non-alcoholic steatohepatitis (NASH), fibrosis, and potentially cirrhosis and hepatocellular carcinoma, posing a serious threat to human health.² Although its pathogenesis involves complex interactions among multiple factors and pathways, dysregulated lipid metabolism, oxidative stress, and inflammatory responses are widely recognized as core drivers of NAFLD initiation and progression.³ Unfortunately, specific therapies capable of effectively reversing NAFLD pathology are currently lacking in the clinic, with interventions primarily focusing on complication management and metabolic parameter regulation. There is an urgent need to dissect the molecular pathogenesis in depth to provide a theoretical basis and novel targets for developing targeted therapies.

Within the complex pathogenic network of NAFLD, mitochondrial dysfunction has garnered increasing attention. As the primary site of cellular energy metabolism and reactive oxygen species (ROS) production, mitochondria play a crucial role in maintaining hepatocyte lipid homeostasis.⁴ Studies indicate that mitochondrial structural and functional damage is closely associated with intrahepatic lipid accumulation, insulin resistance, and inflammation activation. Mitophagy, a highly selective autophagic process responsible for identifying and removing dysfunctional or damaged mitochondria, is a core mechanism for maintaining mitochondrial quality control and cellular homeostasis.⁵ Impaired mitophagy leads to the accumulation of abnormal mitochondria, resulting in excessive ROS generation, energy metabolism disruption, and exacerbated lipotoxicity, thereby promoting hepatocyte steatosis, inflammatory responses, and fibrogenesis.⁶ Consequently, restoring or enhancing mitophagy is considered a promising strategy for halting NAFLD progression,⁷ making the exploration of its regulatory mechanisms highly significant scientifically and clinically.

Although existing research suggests a key role for mitophagy in NAFLD, the related molecular events and their dynamic regulation across different disease stages remain incompletely systematically elucidated. The rapid development of bioinformatics technologies, integrating multi-omics data with artificial intelligence algorithms, offers a novel perspective for mining key disease genes and revealing molecular interaction networks. This study aims to leverage public database resources comprehensively, combining bioinformatics analysis, multiple machine learning models, and reverse network pharmacology approaches to systematically identify mitophagy-related core genes in NAFLD, explore their potential functional mechanisms, and predict corresponding targeted drugs and active compounds. The findings from computational analyses were further validated through in vitro experiments using an NAFLD cell model. The findings are expected to provide new theoretical foundations and data support for understanding disease mechanisms, early diagnosis, and precise treatment of NAFLD.

Furthermore, a critical gap in knowledge lies in the precise mechanism by which defective mitophagy transitions from a cellular metabolic flaw to a driver of systemic inflammation, propelling the progression from benign steatosis (NAFL) to inflammatory steatohepatitis (NASH). It is hypothesized that dysfunctional mitochondria, when not adequately cleared, become a potent source of damage-associated molecular patterns (DAMPs), such as mitochondrial DNA and cardiolipin. These DAMPs activate innate immune sensors (eg., the NLRP3 inflammasome) in hepatocytes and Kupffer cells, triggering the release of pro-inflammatory cytokines (eg., IL-1β, IL-18). This initiates a chemokine gradient that recruits peripheral immune cells, including monocytes and lymphocytes, into the liver parenchyma. Consequently, a self-sustaining cycle of inflammation, cell death, and fibrosis is established. Therefore, identifying the key mitophagy-related genes that malfunction in NAFLD and elucidating their downstream impact on immune cell communication is essential to decode the inflammatory switch in NASH pathogenesis.

Materials and Methods

Data Collection and Preprocessing

Six NAFLD-related datasets were retrieved and downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), comprising 3 microarray datasets, 2 next-generation sequencing datasets, and 10 single-cell RNA sequencing datasets. Details of the datasets are summarized in Table 1. The microarray datasets GSE63067, GSE89632, and GSE164760 were merged to form the training set. Probe annotation and normalization were performed for each dataset individually before merging. The ComBat algorithm was applied to the merged training set data to correct for batch effects, and the effectiveness of batch effect removal was verified by principal component analysis (PCA). The GSE135251 and GSE61260 datasets served as independent validation sets.

Table 1 Overview of Datasets Used in This Study

Mitophagy-Related Gene Set

Genes associated with the biological process “Mitophagy” were retrieved from the Genecards database (https://www.genecards.org/) using the search term “mitophagy”. To ensure high relevance, only genes with a “relevance score” (an internal metric of Genecards reflecting the strength of association with the query term) greater than 1.0 were selected, resulting in a final set of 1765 mitophagy-related genes for subsequent analysis.

Differential Expression Analysis

Differential expression analysis was performed using the R package limma. A linear model was constructed with sample group (healthy vs. NAFLD) as the core variable in the design matrix. Models were fitted using lmFit and eBayes functions with empirical Bayes moderation. Contrasts between groups were calculated using makeContrasts. Significantly differentially expressed genes (DEGs) were defined as those with |log₂ Fold Change| > 0.585 (ie., 1.5-fold change) and an adjusted P-value (FDR) < 0.05.

KEGG Pathway Enrichment Analysis

KEGG pathway enrichment analysis for the identified DEGs was conducted using the clusterProfiler package. Gene symbols were converted to Entrez IDs, and unmatched genes were removed. Enrichment analysis was performed using the enrichKEGG function (parameters: pvalueCutoff=1, qvalueCutoff=1), and significantly enriched pathways were filtered with a threshold of P < 0.05. The top 30 significantly enriched pathways were visualized using bar and dot plots.

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA was used to construct gene co-expression networks. The expression matrix was quality-controlled, filtering out genes with low variance (standard deviation ≤ 0.5). The goodSamplesGenes function checked for and removed samples and genes with excessive missing values. Sample hierarchical clustering was performed to identify and remove outliers. To satisfy the scale-free topology criterion, the pickSoftThreshold function evaluated powers from 1 to 20, selecting the smallest power where the scale-free topology fit index R² first reached above 0.8 as the soft threshold. Using the chosen soft threshold, co-expression modules were constructed via the blockwiseModules function (minModuleSize = 60). Module eigengenes (MEs), defined as the first principal component of each module representing the gene expression profile, were calculated. The correlation between each ME and the binary disease trait (NAFLD = 1, Healthy = 0) was then computed. The module with the highest absolute correlation coefficient (positive or negative) and a statistically significant P-value (< 0.05) was identified as the key module most associated with the disease trait. Module eigengenes were calculated, and highly similar modules (correlation > 0.75) were merged to obtain the final gene module assignments. Key modules most associated with the disease trait (NAFLD vs. healthy) were identified by calculating module-trait relationships.

Machine Learning for Core Gene Screening

Based on the candidate genes from the above analyses, a comprehensive machine learning strategy was employed. A total of 11 mainstream algorithms, including regularized regression, support vector machines, discriminant analysis, and various ensemble learning models, were used. An “feature selection + modeling” two-step process generated 87 model combinations for evaluation. Hyperparameters for all models were optimized via cross-validation. Model performance was comprehensively assessed using accuracy, sensitivity, specificity, F1-score, and the area under the receiver operating characteristic curve (AUC). Results were visualized using multiple methods.

SHAP Model Interpretability Analysis

To interpret the decision-making mechanism of the best-performing model, SHAP (SHapley Additive exPlanations) analysis was performed. An input set based on the identified key features was constructed, and the kernelshap algorithm was used to compute the contribution value (SHAP value) of each feature to the prediction outcome for individual samples. The shapviz package generated various visualizations, including feature importance ranking plots, SHAP summary (beeswarm) plots, SHAP dependence plots, waterfall plots, and force plots for individual samples, revealing feature contribution patterns and interactions from multiple dimensions.

Protein-Protein Interaction (PPI) Network Construction

Core genes identified by machine learning were submitted to the GeneMANIA database (https://genemania.org/) to construct a PPI network, exploring functional associations and potential interactions among the core targets.

Immune Infiltration Analysis

The CIBERSORT algorithm was used to deconvolute transcriptomic data and quantify the relative proportions of 22 immune cell subtypes. The analysis was run with 1000 permutations for significance assessment, and the expression matrix was quantile-normalized to reduce technical bias. Only samples with significant immune cell proportion estimates (P < 0.05) were retained. Differences in immune cell proportions between the NAFLD group and the healthy control group were compared using the Wilcoxon rank-sum test, and results were displayed using stacked bar plots and box plots. Spearman correlation coefficients among immune cell subtypes were calculated and visualized via heatmap. Furthermore, the linkET package was used to analyze associations between core gene expression and immune cell infiltration levels, visualized using correlation network plots.

Gene Set Enrichment Analysis (GSEA)

Single-gene GSEA for the core genes was performed using the clusterProfiler package. Samples were divided into high- and low-expression groups based on the median expression level of the target gene. Genome-wide expression differences between the two groups were calculated and ranked. Custom KEGG gene sets were used for enrichment analysis, significant pathways (P < 0.05) were identified, and enrichment plots were generated.

Single-Cell RNA Sequencing Data Analysis

Data Preprocessing and Quality Control

Single-cell RNA-seq data from 10 samples (5 Healthy, 5 NAFLD) were processed using Seurat (v5). For each sample, expression matrices were read, Seurat objects were created, and quality control thresholds were applied: genes detected in ≥3 cells, cells with detected gene counts between 500 and 7500, unique molecular identifiers (UMIs) > 1000, mitochondrial gene percentage (<20%), and log10(genes/UMI) > 0.8.

Doublet Detection and Filtration

The scDblFinder package was used for doublet identification. Data were normalized and highly variable genes (2000) were selected before converting Seurat objects to SingleCellExperiment objects for doublet scoring. Cells classified as “singlet” with a doublet score < 0.8 were retained.

Batch Effect Correction and Data Integration

The Harmony algorithm was used to correct batch effects across the 10 samples. Data were merged, normalized, and 2000 highly variable genes were selected. Data scaling regressed out the effects of mitochondrial percentage and cell cycle scores (S.Score, G2M.Score). After PCA dimensionality reduction, Harmony integration was run using the first 30 principal components with sample ID as the grouping variable.

Cell Clustering and Annotation

Based on Harmony-corrected dimensions, a K-nearest neighbor graph was constructed using the first 20 dimensions, and cell clustering was performed using the Louvain algorithm (resolution=0.8). Visualization was done via UMAP. Cell types were annotated based on classical marker genes: Monocytes (CD14, LYZ), B cells (CD79A, MS4A1), Endothelial cells (PECAM1, VWF), NK cells (NKG7, GNLY), T cells (CD3D, CD3E), Hepatocytes (ALB, APOA1), Stellate cells (ACTA2, COL1A1), Macrophages (CD68, CD163), etc. Dot plots, violin plots, and feature plots were used to validate annotations.

Pseudotime Trajectory Analysis

Monocle2 was used for trajectory inference. Macrophage/monocyte and T cell subsets were analyzed separately. Seurat objects were converted to cell_data_set objects. After preprocessing (num_dim=50) and UMAP reduction, cells were clustered, and trajectory graphs were learned. Pseudotime values were calculated using the order_cells function, and expression dynamics of key genes (IGF1, MYH11, HYOU1, SPATA18, SCD) along the trajectory were analyzed.

Cell-Cell Communication Analysis

CellChat (v1.6.1) was used to infer cell-cell interactions. Based on the “Secreted Signaling” database, communication probabilities and signaling pathway activity were computed. Circle plots, hierarchy plots, and heatmaps visualized the communication networks and ligand-receptor interactions. Network centrality analysis identified key signaling roles, and non-negative matrix factorization (NMF) was used to identify communication patterns.

Cell Proportion and Differential Expression Analysis

Proportions of different cell types per sample were calculated. The Wilcoxon rank-sum test compared cell proportions between healthy and NAFLD groups. For each cell type, the FindMarkers function was used for intergroup differential expression analysis (parameters: min.pct=0.1, logfc.threshold=0.25) to identify NAFLD-associated, cell type-specific gene expression changes.

Reverse Network Pharmacology and Molecular Docking

Compound-gene interaction data were retrieved from the Comparative Toxicogenomics Database (CTD, https://ctdbase.org/), which was utilized to identify potential pharmaceutical ingredients targeting core genes. The Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) was employed for screening bioactive components, with filtering criteria set as oral bioavailability ≥30% and drug-likeness ≥0.18. Through enrichment analysis of traditional Chinese medicines (TCMs) against the potential pharmaceutical ingredients of core genes, TCMs with statistical significance (P-value < 0.05 and adjusted P-value < 0.05) were selected as candidate herbs targeting the core genes. A “herb-component-target” interaction network was constructed based on the enrichment results, and significant associations were screened using Cytoscape_v3.8.0. Subsequently,

The online molecular docking platform CB-Dock2,⁸ which employs Curvature-based Cavity Detection (CurPocket) for binding site prediction, was used to perform molecular docking between the proteins encoded by the screened core genes and the candidate pharmaceutical ingredients, thereby evaluating their binding affinity and interaction mechanisms.

Experimental Validation of NAFLD Core Genes

To validate the expression differences of NAFLD core genes, HepG2 cells were used as the control group and an in vitro NAFLD model was established. Western Blot and qRT-PCR were used to detect protein and RNA level differences, respectively, and JC-1 staining was used to assess mitochondrial membrane potential changes indicative of mitophagy.

Materials

HepG2 cells were purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences (Shanghai, China). The cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) (Gibco, USA) supplemented with 10% fetal bovine serum (FBS; Invitrogen, USA) and 100 μg/mL penicillin-streptomycin (Suzhou New Cellmei Biotechnology Co., Ltd., China).

The primary antibodies used were as follows: Anti-IGF1 Rabbit Polyclonal Antibody (IGF1, Catalog No. A11985; ABclonal Technology, China), Anti-MYH11 Rabbit Monoclonal Antibody (MYH11, Clone EPR5335; Abcam, Cambridge, UK), Anti-HYOU1 Rabbit Polyclonal Antibody (HYOU1, Catalog No. A1042; ABclonal Technology, Wuhan, China), Anti-SPATA18 Rabbit Polyclonal Antibody (SPATA18, Catalog No. A09906-1; ABclonal Technology, Wuhan, China), Anti-SCD (Stearoyl-CoA Desaturase) Rabbit Monoclonal Antibody (SCD, Catalog No. ab236868; Abcam, Cambridge, UK), and Anti-β-Actin Mouse Monoclonal Antibody (β-Actin, Catalog No. AC026; ABclonal Technology, Wuhan, China).Sodium palmitate (PA) and sodium oleate (OA) were obtained as a ready-made mixture (PA 6 mmol/L, OA 12 mmol/L; Catalog No. KC006) from Xi’an Kunchuang Technology Development Co., Ltd. (China).The mitochondrial membrane potential fluorescent probe JC-1 (Catalog No. C2006) was supplied by Beyotime Biotechnology (China). Carbonyl cyanide m-chlorophenyl hydrazone (CCCP; Product No. C2759) was purchased from Sigma-Aldrich (USA).

Cell Culture

HepG2 cells were cultured in DMEM with 10% FBS and 1% penicillin/streptomycin at 37°C with 5% CO₂.

NAFLD Model Induction

A free fatty acid (FFA) mixture of OA and PA (2:1 molar ratio) complexed with BSA was prepared. HepG2 cells at ~70-80% confluence were treated with 0.5 mM FFA for 24 hours. Control cells received BSA only. Successful lipid accumulation was confirmed by Oil Red O staining.

qRT-PCR

Total RNA was extracted, reverse-transcribed to cDNA, and qPCR was performed using gene-specific primers (Table 2). Relative mRNA levels were calculated using the 2^(-ΔΔCt) method normalized to GAPDH.

Table 2 Primer Sequences for qRT-PCR

Western Blotting

Cells were lysed, proteins quantified (BCA assay), separated by SDS-PAGE, transferred to PVDF membranes, blocked, and incubated with primary antibodies overnight at 4°C, followed by HRP-conjugated secondary antibodies. Signals were detected using ECL chemiluminescence and analyzed with ImageJ software.

JC-1 Staining

Cells were incubated with JC-1 working solution 5 μM at 37°C for 15–20 min, washed, fixed, and imaged by confocal microscopy. JC-1 aggregates (red, high membrane potential) and monomers (green, low membrane potential) were detected. The red/green fluorescence intensity ratio was calculated to assess mitochondrial membrane potential.

All experiments were conducted with a minimum of three independent biological replicates per experimental condition (Control and NAFLD Model). qRT-PCR, JC-1 staining, and Western blot analyses were carried out on these biological replicates. For each biological replicate, three technical replicates were incorporated to ensure data reliability and reduce technical variability.

Statistical Analysis

All bioinformatics analyses were performed using R software (version 4.5.0). Multiple hypothesis testing was corrected using the False Discovery Rate (FDR) method. Group comparisons used appropriate parametric or non-parametric tests. A significance level of P < 0.05 was set. All analyses ensured statistical power through technical and biological replicates.

Results

Data Integration and Mitophagy-Related Gene Expression

Five NAFLD gene expression datasets were obtained from GEO (Table 1). GSE135251 and GSE61260 were designated as validation sets. GSE63067, GSE89632, and GSE164760 were merged to form the training set. Raw data distribution and PCA before normalization showed batch effects (Figure 1A and 1B). After batch effect removal using ComBat and normalization, the integrated training set data showed significantly reduced batch effects (Figure 1C and D).

Figure 1 (A) Boxplot of raw training set data before normalization. (B) PCA plot before batch effect correction. (C) Boxplot of training set data after batch effect correction and normalization. (D) PCA plot after batch effect correction.

Differential Expression and WGCNA Analysis

Differential analysis identified 230 significant DEGs between healthy and NAFLD groups (91 up, 139 down) (Figure 2A). A heatmap displayed the expression patterns of the top 50 DEGs (Figure 2B). KEGG enrichment analysis of DEGs revealed involvement in pathways like Cytokine-cytokine receptor interaction, TNF signaling, Lipid and atherosclerosis, MAPK signaling, and IL-17 signaling (Figure 2C). WGCNA analysis selected a soft threshold power of 10, ensuring a biologically relevant and statistically stable co-expression network (Figure 2D). Co-expression modules were identified; the grey module (320 genes, r = −0.61, p = 8e-18) showed the strongest negative correlation with disease status (Figure 2E). The gene clustering dendrogram and module assignment are shown (Figure 2F). Intersection of the WGCNA key module genes (320) with the mitophagy-related gene set (1765) yielded 17 hub genes (Figure 2G), suggesting a close association of this module with mitochondrial function.

Figure 2 (A) Volcano plot of NAFLD DEGs. (B) Heatmap of top 50 DEGs. (C) Bubble plot of KEGG enrichment for DEGs. (D) Scale-free topology model fit and mean connectivity for different soft-thresholding powers in WGCNA. (E) Module-trait relationships heatmap. (F) Gene clustering dendrogram and module colors. (G) Venn diagram showing overlap between WGCNA grey module genes and mitophagy-related genes (17 hub genes).

Machine Learning Identifies Core Genes

Based on the 17 hub genes, multiple machine learning models were evaluated across three cohorts (training set, GSE135251, GSE61260) (Figure 3A). The glmBoost+GBM combination performed best, identifying 8 key genes: IGF1, MYH11, HYOU1, SPATA18, EPS8L2, SCD, ACACA, FEN1. Among these, IGF1, MYH11, HYOU1, SPATA18, and SCD were also significant DEGs, as marked on the volcano plot (Figure 3B). This model showed excellent diagnostic performance with AUCs of 0.974, 0.951, and 0.756 in the training, GSE135251, and GSE61260 cohorts, respectively. Confusion matrices demonstrated robust performance across cohorts (Figure 3C–E). Analysis of feature selection frequency highlighted MYH11, SCD, and SPATA18 as the most consistently selected genes, indicating their importance in disease discrimination (Figure 3F).

Figure 3 (A) AUC performance of various machine learning models across cohorts. (B) Volcano plot highlighting the 8 key genes identified by glmBoost+GBM; IGF1, MYH11, HYOU1, SPATA18, SCD are significant DEGs. (C–E) Confusion matrices for the training set (C), GSE135251 validation set (D), and GSE61260 validation set (E). (F) Bar plot showing the selection frequency of the most stable features during model building.

SHAP Analysis for Model Interpretation

SHAP analysis revealed that high SHAP values for IGF1 were concentrated, indicating its significant positive influence on model output (predicting NAFLD) (Figure 4A). SPATA18 and HYOU1 showed moderate positive influences, while ACACA and FBN1 had relatively smaller, mainly negative contributions. The mean absolute SHAP value ranking confirmed IGF1 (0.665) as the most important feature, followed by HYOU1 (0.450) and SPATA18 (0.334) (Figure 4B). For an individual NAFLD sample prediction (output value 2.26), IGF1 contributed +0.699, being the most important diagnostic factor (Figure 4C). Conversely, for a healthy sample prediction (output value −0.463), IGF1 was the primary negative contributor (−1.21) (Figure 4D). Bar plots for specific samples further emphasized IGF1’s dominant role (Figure 4E and F). SHAP dependence plots suggested a positive correlation between HYOU1 and IGF1 SHAP values, and a potential interaction between MYH11 and IGF1. An inverse relationship was observed between SPATA18 and SCD (Figure 4G).

Figure 4 SHAP analysis. (A) Bee swarm plot of SHAP values for main features. (B) Feature importance based on mean absolute SHAP values. (C and D) SHAP force plots decomposing predictions for an individual NAFLD sample (C) and a healthy sample (D). (E and F) Bar plots showing feature contributions for specific NAFLD (E) and healthy (F) samples. (G) SHAP dependence plots for key feature pairs.

Core Gene Relationships and GSEA

The five final core genes—IGF1, MYH11, HYOU1, SPATA18, and SCD—showed significant expression differences between disease and healthy groups. IGF1, MYH11, and HYOU1 were downregulated, while SPATA18 and SCD were upregulated in NAFLD (Figure 5A). Correlation analysis revealed a significant negative correlation between IGF1 and SCD (Figure 5B). PPI network analysis positioned IGF1 centrally, connecting it to processes like IGF-like growth factor signaling, oxygen metabolism, and T-cell proliferation (Figure 5C). GSEA for the five core genes indicated their involvement in metabolic regulation, MAPK signaling, endoplasmic reticulum stress, and fibrosis pathways (Figure 5D).

Figure 5 (A) Expression differences of core genes between healthy and NAFLD groups (**P < 0.01, ***P < 0.001). (B) Correlation analysis among core genes (*P < 0.05,**P < 0.01, ***P < 0.001). (C) PPI network and functional associations of core genes. (D) GSEA enrichment plots for core genes.

Immune Infiltration Landscape

Stacked bar plots and box plots showed distinct immune cell composition between groups. The NAFLD group had relatively higher proportions of T cells gamma delta, Monocytes, Macrophages M2, Dendritic cells resting, and Mast cells resting. Notably, T cells CD4 memory activated and Mast cells activated were significantly increased in NAFLD (Figure 6A and B). Correlation heatmap revealed significant positive and negative correlations among various immune cell types (Figure 6C). Core gene expression was correlated with immune cell infiltration levels. For instance, IGF1 negatively correlated with T cells CD4 naive but positively with T cells gamma delta. MYH11 negatively correlated with Macrophages M1. HYOU1 showed negative correlations with T cells CD4 memory activated, T cells gamma delta, and Macrophages M2. SCD negatively correlated with B cells naive and Dendritic cells activated but positively with T cells CD4 naive and Macrophages M2 (Figure 6D).

Figure 6 Immune infiltration analysis. (A) Stacked bar plot of immune cell proportions. (B) Box plots comparing immune cell proportions between groups (*P < 0.05, **P < 0.01, ***P < 0.001). (C) Heatmap of correlations among immune cell types. (D) Network plot showing correlations between core gene expression and immune cell infiltration levels.

Single-Cell Transcriptomics Reveals Cellular Heterogeneity

UMAP visualization identified 20 distinct cell clusters from the integrated single-cell data (Figure 7A). These were annotated into major liver cell types: T cells, Monocytes, NK cells, Endothelial cells, Macrophages, B cells, and Stellate cells (Figure 7B). Proportion analysis showed a significant increase in T cells and Monocytes, and a decrease in Endothelial cells and B cells in NAFLD (Figure 7C and E). Dot plot of marker gene expression supported the annotation (Figure 7D). Differential expression analysis within each cell type identified numerous NAFLD-associated genes (Figure 7F). Feature plots showed specific expression patterns of the core genes: IGF1 and SCD were concentrated in Macrophages, MYH11 in Stellate cells, SPATA18 in Hepatocytes, and HYOU1 was broadly expressed (Figure 7G).

Figure 7 Single-cell analysis. (A) UMAP of 20 cell clusters. (B) UMAP annotated by major cell types. (C) Bar plot of cell type proportions per sample/group. (D) Dot plot of marker gene expression for cell type annotation. (E) Box plots comparing cell type proportions between groups. (F) Volcano plots of DEGs for 8 major cell types. (G) UMAP feature plots showing expression distribution of core genes.

Pseudotime Trajectory Analysis

Monocytes/Macrophages: Trajectory analysis revealed three branches, indicating heterogeneity (Figure 8A and B). Pseudotime analysis showed dynamic changes along the trajectory (Figure 8C). IGF1 and SCD were upregulated during the monocyte-to-macrophage transition (Figure 8D).

Figure 8 Pseudotime analysis. (A–D) Monocytes/Macrophages: (A) UMAP, (B) State classification, (C) Pseudotime trajectory, (D) Marker gene expression. (E–H) Hepatocytes: (E) UMAP by group, (F) State classification, (G) Pseudotime trajectory, (H) Marker gene expression. (I–L) T cells: (I) UMAP by group, (J) State classification, (K) Pseudotime trajectory, (L) Marker gene expression.

Hepatocytes: Distinct spatial distribution was observed between healthy and NAFLD hepatocytes (Figure 8E). Three hepatocyte subpopulations were identified (Figure 8F). The trajectory depicted a transition from healthy to pathological states (Figure 8G). Core genes showed state-specific expression changes (Figure 8H).

T cells: Distribution differences were noted between healthy and NAFLD T cells (Figure 8I). Multiple T cell states were identified (Figure 8J). The trajectory revealed T cell evolution (Figure 8K). MYH11, HYOU1, and SCD were notably upregulated in activated states (Figure 8L).

Cell-Cell Communication Analysis

The circle plot showed extensive interactions, particularly involving Macrophages, Monocytes, and T cells (Figure 9A). Heatmaps indicated that T cells and Macrophages had strong outgoing and incoming signaling for immune-related pathways (eg., CCL, MIF, VISFATIN). Hepatocytes showed strong outgoing signals for MIF, SPP1 and incoming signals for MK, PTN (Figure 9B). River plots showed specific signal pattern distributions across cell types (Figure 9C and D). The MIF signaling pathway involved extensive communication from Hepatocytes to immune cells (Figure 9E and F). The MK pathway showed signaling from Endothelial cells, Stellate cells, and Hepatocytes to Hepatocytes (Figure 9G and H).

Figure 9 Cell-cell communication. (A) Circle plots showing number (top) and strength (bottom) of interactions. (B) Heatmaps of outgoing (left) and incoming (right) signal patterns per cell type. (C and D) River plots for outgoing (C) and incoming (D) signals. (E and F) MIF signaling: chord diagram (E) and heatmap of ligand-receptor pair strengths (F). (G and H) MK signaling: chord diagram (G) and heatmap of ligand-receptor pair strengths (H).

Reverse Network Pharmacology and Molecular Docking

Drug enrichment analysis constructed a network linking 14 significant herbs, 3 active components, and 2 core targets (IGF1, SCD). Twelve herbs acted on IGF1 via baicalein, Niuhuang acted on IGF1 and SCD via Deoxycholic Acid, and Bajiaolian acted on SCD via deoxypodophyllotoxin (Figure 10A). Molecular docking showed strong binding potential: IGF1-Baicalein (Vina score: −6.1), IGF1-Deoxycholic Acid (−6.4), SCD-Deoxycholic Acid (−7.3), SCD-Deoxypodophyllotoxin (−7.3) (Figure 10B–E).

Figure 10 (A) Network of herb-component-target interactions. (B–E) Molecular docking results for (B) IGF1 and Baicalein, (C) IGF1 and Deoxycholic Acid, (D) SCD and Deoxycholic Acid, (E) SCD and Deoxypodophyllotoxin.

Experimental Validation

At the transcriptional level, quantitative real-time polymerase chain reaction (qRT-PCR) analysis (Figure 10B) validated the significant dysregulation of the identified core genes. Consistent with our bioinformatic predictions, IGF1, MYH11, and HYOU1 expression was markedly downregulated, whereas SPATA18 and SCD were upregulated in the FFA-induced NAFLD cell model compared with control groups (P < 0.05 for all comparisons; Figure 11A). These changes were further validated at the protein level through Western Blot analysis, which revealed consistent expression patterns of the core proteins (Figure 11B). The Western Blot bands confirmed that the protein expression levels exhibited concordant trends with the mRNA expression profiles observed in qRT-PCR (Figure 11A), thereby providing robust evidence for the translational regulation of the identified core genes. Oil Red O staining (Figure 11C) revealed that, compared to the control group, cells treated with FFA exhibited significant red staining, indicating intracellular lipid accumulation. JC-1 staining showed high red/green fluorescence ratio (high membrane potential, low mitophagy) in control cells, an increased green fluorescence (decreased membrane potential, increased mitophagy) in some NAFLD cells, and predominant green fluorescence (low membrane potential, high mitophagy) in the CCCP positive control (Figure 11D).

Figure 11 Experimental validation. (A) qRT-PCR results for core gene mRNA expression (**P < 0.01, ***P < 0.001). (B) Western Blot bands for core proteins. (C) Oil Red O staining of control and NAFLD model HepG2 cells. (D) JC-1 fluorescence staining (red: aggregates, high potential; green: monomers, low potential).

Discussion

NAFLD is a chronic liver disease closely associated with metabolic disorders. Its global prevalence is continuously rising, making it a significant public health challenge. The disease spectrum progresses from simple hepatic steatosis to non-alcoholic steatohepatitis (NASH), liver fibrosis, and even cirrhosis and hepatocellular carcinoma, posing a serious threat to human health. Although dysregulated lipid metabolism, oxidative stress, and inflammatory responses are considered key links in the pathogenesis of NAFLD, a systematic analysis of its molecular mechanisms is insufficient, and the regulatory role of mitophagy in the pathological process has not been fully elucidated. This study aims to provide new perspectives for investigating disease mechanisms and developing targeted intervention strategies by integrating multi-omics data, applying machine learning algorithms, and reverse network pharmacology to systematically screen for mitophagy-related core genes in NAFLD.

Through bioinformatic analysis of multiple GEO cohorts, we first identified 230 genes significantly differentially expressed between healthy and NAFLD groups. Further combined with Weighted Gene Co-expression Network Analysis (WGCNA), we screened a gene module significantly associated with disease status, which overlapped with the mitophagy gene set, yielding 17 overlapping genes. Building on this, we employed a comprehensive strategy involving 11 machine learning algorithms for feature selection and ultimately identified five core genes: IGF1, MYH11, HYOU1, SPATA18, and SCD. These genes demonstrated excellent disease discrimination performance across multiple validation cohorts (maximum AUC reaching 0.974). Functional enrichment analysis revealed that these genes were significantly enriched in key biological processes such as the PI3K-Akt signaling pathway, endoplasmic reticulum stress response, mitophagy, and lipid metabolism, suggesting their potential involvement in the pathogenesis of NAFLD by coordinating metabolic homeostasis and cellular quality control mechanisms.

IGF1 is a polypeptide hormone produced primarily by the liver in response to growth hormone (GH) stimulation. It is widely expressed throughout the body and is a core member of the insulin-like growth factor system, exerting potent mitogenic, anti-apoptotic, and metabolic regulatory effects by activating its receptor IGF1R (structurally similar to the insulin receptor).⁹ Numerous clinical and basic studies¹⁰ have shown that circulating and local hepatic IGF1 levels are significantly reduced in NAFLD patients, and the degree of reduction is negatively correlated with the severity of NAFLD, consistent with our finding of significant downregulation of IGF1 in the NAFLD group. The decrease in IGF1 not only aggravates insulin resistance but also promotes impaired fatty acid β-oxidation, leading to lipid droplet accumulation in hepatocytes, forming the molecular basis for NAFLD development.¹¹ Furthermore, the availability of dietary carbohydrate-to-protein calorie ratio affects hepatic fat deposition through the IGF-1/PI3K/Akt signaling pathway.¹² Excessive accumulation of free fatty acids, sphingolipids, and lipid peroxidation (lipotoxicity) leads to hepatocyte injury, resulting in the production of damage-associated molecular patterns (DAMPs), which activate hepatic sinusoidal endothelial cells and resident macrophages (Kupffer cells) to promote inflammation and attract immune cells. Additionally, DAMPs and inflammatory cytokines activate hepatic stellate cells (HSCs), leading to remodeling of the local extracellular matrix and deposition of excess collagen and other proteins, causing fibrosis. IGF1, by activating the PI3K-Akt-mTOR signaling pathway, is an important sensor of cellular energy and nutrient status.¹³ It is well known that mTOR is a key negative regulator of autophagy, including mitophagy.¹⁴ Under nutrient-sufficient conditions, the active IGF1-Akt-mTOR pathway suppresses autophagy to facilitate anabolism. However, in the NAFLD state, decreased IGF1 levels may lead to aberrant mTOR activity, which disrupts normal metabolic homeostasis and affects the baseline level of mitophagy.¹⁵ Furthermore, IGF-1 prevents the accumulation of autophagic vesicles and increased cell death, as well as promotes the fusion and degradation rates of autophagosomes to lysosomes, thereby contributing to cell survival.¹⁶ IGF1 deficiency leads to insulin resistance and lipid overload, subjecting mitochondria to immense β-oxidation pressure and generating excessive ROS. ROS is a major inducer of mitochondrial damage, and mitochondrial damage is a key signal for initiating mitophagy.¹⁷ Therefore, IGF1 downregulation, excessive mitochondrial ROS production, and mitochondrial damage should theoretically trigger mitophagy. However, under persistent metabolic stress, this compensatory mitophagy may be insufficient or dysfunctional, failing to effectively clear damaged mitochondria, leading to the accumulation of “zombie” mitochondria, further amplifying oxidative stress and cellular injury, forming a vicious cycle.

MYH11 encodes smooth muscle myosin heavy chain, a core component of the contractile apparatus in vascular and smooth muscle cells.¹⁸ It is a marker gene for mature vascular smooth muscle cells (VSMCs), and its expression is crucial for maintaining vascular tension and structure.¹⁹ In the liver, MYH11 expression is primarily associated with the vascular system and activated hepatic stellate cells (HSCs).²⁰ During NAFLD/NASH progression, pathological changes occur in the liver microvasculature, a process termed “sinusoidal capillarization”,²¹ where hepatic sinusoidal endothelial cells lose their fenestrae, form a basement membrane, and cells expressing smooth muscle markers appear. The downregulation of MYH11 may reflect dysfunction and abnormal remodeling of the liver vascular system, affecting hepatic blood supply and oxygenation, and exacerbating hepatocyte hypoxia and injury.²² It is noteworthy, however, that quiescent HSCs do not express MYH11. When the liver is injured, HSCs are activated and transform into myofibroblasts, acquire contractility, begin to express α-SMA (smooth muscle actin) and MYH11, and secrete large amounts of collagen, leading to liver fibrosis.²³ Therefore, the detection of MYH11 downregulation in whole liver tissue might imply: 1) The detected signal primarily originates from hepatocytes or vascular cells, rather than activated HSCs; or 2) In late-stage disease, the activation state of HSCs might undergo further changes. The altered expression of MYH11 is more likely a downstream marker of NAFLD-associated vascular pathology and fibrotic progression. Mitophagy dysregulation occurs within hepatocytes as a result of metabolic injury; the downregulation of MYH11 indicates a deteriorating liver microenvironment, indirectly reflecting the stress level endured by mitochondria.

HYOU1, also known as ORP150 or GRP170, is an important molecular chaperone in the endoplasmic reticulum (ER) lumen, belonging to the HSP70 family.²⁴ It plays a key role in protein folding, assembly, and transport, and is strongly induced under stress conditions such as hypoxia and nutrient deprivation, serving as a core response protein for endoplasmic reticulum stress (ERS).²⁵ Persistent ERS is a hallmark of NAFLD. Lipid overload (particularly saturated fatty acids) and oxidative stress disrupt ER homeostasis.²⁶ During the initial stages of stress, HYOU1, as a protective protein, is typically upregulated to assist in the refolding of misfolded proteins, alleviating ERS, and maintaining cell survival.²⁷ The observed downregulation of HYOU1 in the NAFLD group is a critical signal! It indicates that persistent and severe ERS has overwhelmed the cell’s adaptive response mechanisms, and the compensatory capacity of the chaperone system is exhausted. This exhaustion signifies that the unfolded protein response (UPR) shifts from a pro-survival mode to a pro-apoptotic mode, significantly increasing the risk of hepatocyte death.²⁸ The ER and mitochondria are closely connected via mitochondria-associated ER membranes (MAMs). Persistent ERS leads to calcium ion (Ca²⁺) leakage from the ER. Excessive mitochondrial Ca²⁺ uptake causes membrane potential collapse and ROS bursts, severely damaging mitochondria.²⁹ As a key chaperone protein in the MAM region, the downregulation of HYOU1 indicates severe ER insufficiency, unable to effectively manage protein folding and Ca²⁺ homeostasis.³⁰ This exacerbates the mitochondrial damage initiated by ERS. The presence of numerous damaged mitochondria potently activates mitophagy pathways (eg., PINK1-Parkin) in an attempt to clear them.³¹ Therefore, HYOU1 downregulation does not inhibit autophagy but rather, by exacerbating mitochondrial damage, indirectly increases the burden and necessity of mitophagy. It is an early warning signal of cellular collapse.

SPATA18, also known as Mieap (Mitochondria-eating protein), is a core gene that directly executes mitophagy. Its encoded protein is recruited to the mitochondrial outer membrane in response to mitochondrial damage, promoting the recognition and engulfment of damaged mitochondria by lysosomes.³² In the NAFLD environment, due to factors such as lipotoxicity, ROS, and ERS (as indicated by HYOU1 downregulation), mitochondria suffer extensive damage.³³ The upregulation of SPATA18 expression represents an active, compensatory cellular response. The cell attempts to enhance its mitophagic capacity to accelerate the clearance of these dysfunctional mitochondria, thereby maintaining mitochondrial network quality and cell survival. SPATA18/Mieap is a key effector protein in the p53-dependent mitophagy pathway, and its upregulation directly demonstrates the activation of the mitophagy pathway in the NAFLD liver.³⁴ The expression level of SPATA18 can, in a sense, reflect the severity of mitochondrial damage within the cell. Its significant upregulation indicates that the mitochondrial quality control system is under tremendous pressure, engaged in a “clearance battle.”

SCD is a key enzyme in lipid metabolism that catalyzes the desaturation of saturated fatty acids (SFAs) (eg., palmitic acid, stearic acid) to generate monounsaturated fatty acids (MUFAs) (eg., palmitoleic acid, oleic acid). This is a necessary step for the synthesis of triglycerides (TGs), phospholipids, and cholesteryl esters.³⁵ Saturated fatty acids (SFAs) possess strong lipotoxicity and can induce hepatocyte death by activating mitochondrial death pathways and increasing ER stress.³⁶ SCD, by converting SFAs to MUFAs, reduces the cellular lipotoxicity level,³⁷ representing a self-protective mechanism. However, MUFAs are the preferred substrates for triglyceride (TG) synthesis.³⁸ Excessively high SCD activity provides ample raw materials for TG synthesis, conversely promoting lipid accumulation within hepatocytes, ie., exacerbating hepatic steatosis. Therefore, the upregulation of SCD is a “double-edged sword.” By reducing SFA accumulation, SCD upregulation alleviates mitochondrial damage caused by lipotoxicity. SFAs can directly damage mitochondrial membranes and function, whereas MUFAs are less toxic. From this perspective, SCD upregulation, by improving the lipid environment, indirectly reduces the trigger signals for mitophagy. However, the MUFAs produced by SCD can be incorporated into the phospholipids of the mitochondrial membrane itself, altering its fluidity, thereby affecting mitochondrial function (eg., electron transport chain efficiency) and the efficiency of membrane fusion events during mitophagy.³⁹ The upregulation of SCD represents the cell’s effort to cope with lipotoxic stress. By converting harmful SFAs into relatively harmless MUFAs and storing them, it may temporarily alleviate the assault on mitochondria in the short term, but in the long run, it “adds fuel to the fire” of hepatocyte steatosis. Its relationship with mitophagy is one of indirect regulation, primarily achieved by altering the level of lipid stress experienced by the mitochondria.

This study employed an integrative multi-omics and machine learning approach to identify and validate five core genes—IGF1, MYH11, HYOU1, SPATA18, and SCD—intricately linked to mitophagy and NAFLD pathogenesis. Their expression patterns paint a coherent picture of the disease’s molecular landscape. Downregulation of IGF1, a key metabolic and anti-inflammatory hormone, sets the stage for insulin resistance and reduced antioxidant capacity, priming hepatocytes for lipid accumulation and oxidative stress. This metabolic dysregulation, particularly lipotoxicity from saturated fatty acids, induces endoplasmic reticulum stress. The observed downregulation of HYOU1, an ER chaperone, signals the exhaustion of the protective unfolded protein response, pushing cells towards apoptosis and exacerbating mitochondrial damage via ER-mitochondria coupling. This severe mitochondrial injury triggers a compensatory response, evidenced by the upregulation of SPATA18 (Mieap), a direct executor of mitophagy, in an attempt to clear dysfunctional organelles. Concurrently, the upregulation of SCD represents a double-edged sword; it mitigates direct lipotoxicity by converting SFAs to MUFAs but also fuels hepatic steatosis by providing substrates for triglyceride synthesis. Finally, the downregulation of MYH11, a smooth muscle marker, likely reflects pathological vascular remodeling (“sinusoidal capillarization”) and alterations in hepatic stellate cell phenotype, contributing to the progression of fibrosis and microenvironment disruption. This cascade ultimately creates a pro-inflammatory, pro-fibrotic milieu, as reflected in our immune infiltration and single-cell analyses.

Our data reveal a characteristic NAFLD immune signature: a state of chronic, low-grade inflammation and aberrant immune regulation. This is not a robust acute inflammatory response but rather a state of long-term chronic immune dysregulation. The elevated proportions of T cells and decreased proportions of Macrophages in the disease group suggest that the immune system is skewed towards a suppressive, tolerant, or chronic response pattern, which favors disease persistence rather than clearance.

IGF1 is not only a metabolic hormone but also a significant immunomodulatory factor.⁴⁰ It promotes the proliferation, differentiation, and survival of various immune cells (T cells, B cells, macrophages).⁴¹ Importantly, IGF1 signaling possesses anti-inflammatory properties, inhibiting macrophage polarization towards the pro-inflammatory M1 phenotype and enhancing the function of regulatory T cells (Tregs). The deficiency of IGF1 disrupts this immune balance. The attenuation of its anti-inflammatory effect creates conditions conducive to M1 macrophage activation. Concurrently, its inability to adequately support healthy T cell responses may lead to inefficient immune reactions, aligning with the observed state of adaptive immunity in NAFLD.⁴² The IGF1 signaling pathway is well-characterized for its role in inhibiting macrophage polarization toward the pro-inflammatory M1 phenotype and enhancing the immunosuppressive functions of regulatory T cells.⁴³ Consequently, downregulation of this pathway may directly attenuate the suppression of M1 polarization and impair immunoregulatory mechanisms, collectively fostering a pro-inflammatory microenvironment.

Although MYH11 itself is not an immune molecule, its downregulation, as a marker of vascular smooth muscle, indicates pathological remodeling of the liver’s vascular system (sinusoids), specifically “sinusoidal capillarization.” This alters the hepatic microcirculation, affecting the infiltration and trafficking of immune cells. Sinusoidal structural abnormalities can lead to hypoxia and accumulation of cellular debris; these “danger signals” can, in fact, recruit and activate pro-inflammatory M1 macrophages.⁴⁴ Functional vasculature is essential for dendritic cell (DC) migration. Disordered vascular structure may impede DC activation and migration (manifested as reduced activated DCs), weakening antigen presentation and the initiation of adaptive immunity.⁴⁵ Whole-organ transcriptome sequencing captures the average transcriptional profile across all constituent cells. The observed downregulation of MYH11 may suggest that vascular-derived signals—originating from portal or central venous smooth muscle cells—are diluted by transcriptional contributions from highly proliferating hepatocytes, immune cells, or steatotic hepatocytes. MYH11 expression may exhibit dynamic variation across different stages of fibrotic activation. The study cohort likely includes a substantial proportion of patients with early-stage steatosis, in which hepatic stellate cells (HSCs) may not yet be fully activated into the high-MYH11-expressing myofibroblast state.⁴⁶ Alternatively, in advanced disease, persistent inflammatory milieus could lead to dysregulation of vascular smooth muscle cells and loss of characteristic markers.⁴⁷ Activated HSCs (myofibroblasts) may undergo further differentiation into alternative states, contributing to the downregulation of MYH11. This pattern may reflect a phenotypic transition from a contractile (high MYH11) to a pro-inflammatory and pro-fibrotic (potentially low MYH11) state.⁴⁸

HYOU1 is a guardian against endoplasmic reticulum stress (ERS). ERS is a potent immune modulator. Unresolved ERS can induce apoptosis and inflammatory responses through UPR signaling.⁴⁹ Severe ERS can directly suppress T cell activation and function, even inducing T cell apoptosis,⁵⁰ which explains the weakened effector T cell response.

SPATA18-mediated mitophagy is a process for clearing intracellular “debris.” However, when this process becomes overactivated, it can be associated with immunogenic cell death (ICD).⁵¹ ICD releases damage-associated molecular patterns (DAMPs), such as ATP and mitochondrial DNA, which can activate the innate immune system and recruit antigen-presenting cells. The increase in SPATA18 is a direct response to severe mitochondrial damage. It may indirectly contribute to the infiltration and activation of innate immune cells, such as the observed elevation in Monocytes and M2 macrophages, by promoting ICD and DAMP release. These cells are recruited to clear dead cell debris, but their persistent presence leads to chronic inflammation and fibrosis.⁵²

SCD alters the lipid composition, and specific lipids are potent immune signaling molecules.⁵³ SCD converts saturated fatty acids (SFAs) into monounsaturated fatty acids (MUFAs). SFAs (eg., palmitic acid) can activate the NLRP3 inflammasome and promote M1 polarization,⁵⁴ whereas MUFAs are considered more likely to promote M2 polarization and are less inflammatory. SCD is a central executor of lipid-mediated immune regulation. Its upregulation directly shapes an anti-inflammatory/pro-fibrotic (M2-type) immune microenvironment by altering lipid metabolites, which is a core mechanism for NAFLD persistence. From a therapeutic perspective, these genes provide multiple targets: supplementing IGF1 or enhancing its signaling could simultaneously improve metabolic and immune status. Inhibiting SCD might disrupt pro-fibrotic M2 polarization, but caution is needed to avoid exacerbating lipotoxicity. Alleviating ERS could reduce cellular stress and the consequent inflammatory signaling. Briefly, the upregulation of SCD alters the lipid composition by elevating the proportion of monounsaturated fatty acids, exemplified by oleic acid.⁵⁵ These lipid metabolites may serve as ligands for peroxisome proliferator-activated receptors (PPARs).⁵⁶ Specifically, the PPAR-γ signaling pathway acts as a critical regulator driving macrophage polarization toward the anti-inflammatory M2 phenotype.⁵⁷ This establishes a plausible mechanistic connection to the observed increase in M2 macrophages.

Traditional bioinformatics analyses often focus on single-factor exploration and struggle to integrate multi-dimensional omics data. Machine learning models, particularly multi-model ensemble and cross-validation techniques, can effectively enhance the accuracy of feature selection and model generalizability. This study, by combining 11 mainstream machine learning methods through a “feature selection + model fusion” pipeline, obtained a stable combination with excellent diagnostic performance. SHAP interpretability analysis not only revealed the contribution of individual core genes to the model but also intuitively demonstrated the mechanism by which multiple features synergistically influence disease discrimination. This approach breaks through the limitations of single biomarkers, enabling multi-level and multi-angle disease prediction, and provides a new tool for clinical risk stratification and individualized intervention.

The application of multiple machine learning algorithms provided a robust framework for feature selection, minimizing bias and identifying the most stable and informative gene signatures. The high AUC values across independent validation cohorts underscore the diagnostic potential of this gene panel. The SHAP analysis further enhanced the interpretability of the model, clearly delineating the individual contributions of each gene and revealing potential synergistic or antagonistic interactions, such as the inverse relationship between SPATA18 and SCD.

The single-cell and cell communication analyses provided unprecedented resolution, localizing core gene expression to specific cell types (eg., SPATA18 in hepatocytes, IGF1/SCD in macrophages) and revealing dynamic changes along cellular trajectories. This highlights the cell-type-specific roles of these genes and the complex intercellular crosstalk, such as hepatocyte-derived MIF signaling, that orchestrates the NAFLD microenvironment.

Finally, the reverse network pharmacology and molecular docking suggest promising repurposing opportunities for traditional medicine. The interactions of baicalein, deoxycholic acid, and deoxypodophyllotoxin with IGF1 and SCD support investigating them as potential therapies for identified pathways. Computational docking predicted favorable binding and affinity, indicating promise for drug development. However, these in silico projections require experimental validation via in vitro assays and functional analyses to confirm efficacy.

Limitations and Future Directions

This study is subject to several limitations. First, the in vitro validation was conducted using the HepG2 hepatoma cell line. Although HepG2 cells are widely employed in modeling NAFLD, they exhibit inherent metabolic disparities compared to primary human hepatocytes, which may affect the biological relevance and translational potential of the observed molecular alterations. Future investigations should utilize primary hepatocytes or more physiologically representative in vitro models, such as hepatocyte spheroids or co-culture systems, to reinforce these findings. Additionally, the reliance on publicly available databases introduces potential batch effects and population heterogeneity; thus, prospective validation utilizing multi-center cohorts is required. The correlations identified between core genes and immune infiltration are derived solely from transcriptomic data and necessitate confirmation at the protein level, complemented by functional assays. Moreover, the predicted drug–target interactions require experimental validation through both in vitro and in vivo studies. Although the machine learning models exhibited strong predictive accuracy, their generalizability must be assessed in larger and more diverse populations. Finally, the intricate mechanisms through which these core genes regulate mitophagy merit further detailed molecular investigation.

Future work should focus on: 1) Collecting large-scale, multi-omics clinical cohorts; 2) Integrating proteomic and metabolomic data for a systems-level understanding; 3) Employing spatial transcriptomics to resolve cellular interactions; 4) Conducting high-throughput drug screens and functional validation of predicted compounds; 5) Developing clinically applicable diagnostic tools based on the core gene panel; and 6) Elucidating the precise molecular mechanisms by which these core genes regulate mitophagy and immune cell function.

Conclusion

In conclusion, this multi-omics study identifies IGF1, MYH11, HYOU1, SPATA18, and SCD as central regulators linking mitophagy deficiency to the inflammatory immune microenvironment in NAFLD. Our work not only provides a novel diagnostic gene signature but also illuminates specific molecular pathways and candidate compounds that could be targeted to disrupt the vicious cycle of metabolic stress, inflammation, and fibrosis in NASH.