Introduction
Liver fibrosis is a pathological condition present in most chronic liver diseases, characterized by excessive proliferation and deposition of extracellular matrix (ECM) within liver tissue. This process leads to abnormal alterations in liver architecture and impairs normal liver function. Fundamentally, liver fibrosis represents an excessive repair response following chronic liver injury, which can become irreversible over time. Persistent fibrosis, accompanied by the necrosis and apoptosis of normal hepatocytes, results in the accumulation of ECM. Eventually, this ECM forms scar tissue, progressively replacing the liver parenchyma and leading to cirrhosis, portal hypertension, or even liver cancer, ultimately resulting in liver failure.1 Chronic liver diseases, including those caused by hepatitis viruses, alcohol, drugs, toxins, parasites, metabolic and genetic factors, cholestasis, and immune abnormalities, are defined by a disease course lasting more than six months. Long-term pathogenic stimulation, abnormal metabolism, and immune-inflammatory responses can damage hepatocytes and trigger the onset of liver fibrosis. Therefore, liver fibrosis is a common outcome in most chronic liver diseases,2 and its further progression can lead to cirrhosis, severely impacting health and life expectancy. Although reliable data on the incidence and prevalence of liver fibrosis are lacking, epidemiological studies and clinical pathology suggest that liver fibrosis is prevalent globally. More than 1.5 billion people worldwide suffer from chronic liver diseases, with the most common causes being non-alcoholic fatty liver disease, hepatitis B virus (HBV), hepatitis C virus (HCV), and alcoholic liver disease.3 Cirrhosis, a direct consequence of liver fibrosis, accounts for 1.82% of the global disease burden, resulting in 1.2 million deaths annually.4 Studies have shown that etiological treatments in Western medicine, such as long-term antiviral therapy for HBV, can effectively suppress viral replication and reverse chronic hepatitis B-related liver fibrosis.5 However, current treatments targeting the underlying etiology have limitations in fully suppressing inflammation, and once the mechanism of fibrosis is triggered, it tends to progress actively. Therefore, antifibrotic therapies targeting the proliferation and degradation of fibrous tissue are crucial and constitute an important therapeutic approach for chronic liver disease. In chronic liver diseases lacking effective etiological treatments, antifibrotic therapy becomes even more critical. Cirrhosis, the end stage of liver fibrosis, underscores the importance of preventing or reversing liver fibrosis as one of its foundational therapeutic strategies.6
Due to the pathological mechanisms of liver fibrosis, which involve a systemic reparative response following liver damage, drugs targeting single molecules often fail to show clinical efficacy. Consequently, there are currently no clearly effective chemical or biological drugs available for clinical use. Traditional Chinese medicine (TCM) is a comprehensive medical system with a robust theoretical framework derived from millennia of clinical practice. Recent decades of research and clinical application have demonstrated that TCM is effective in treating liver diseases, particularly in the prevention and management of liver fibrosis.7 Although TCM did not originally recognize liver fibrosis as a distinct clinical entity, the pathological changes and clinical manifestations of liver fibrosis (including cirrhosis) have long been classified under terms such as “accumulation” or “pain in the hypochondrium.” This understanding has been widely accepted in clinical practice for decades. The primary etiology of liver fibrosis varies, and clinical manifestations differ; however, the fundamental pathogenesis is often one of deficiency and stagnation, with long-term retention of toxins that impair liver meridians, causing qi stagnation and blood stasis, which can be summarized as “deficiency leading to accumulation.” Liver fibrosis essentially results from damage to liver parenchyma and deficiency of the liver and kidney yin, leading to an inability to transform qi, resulting in impaired circulation of qi and blood, blood stasis, and the formation of accumulations.8 The “deficiency” primarily manifests as spleen qi deficiency, liver qi deficiency, and liver-kidney yin deficiency. Qi deficiency reflects functional damage and decline, while liver-kidney yin deficiency denotes deeper pathological changes in liver parenchyma. Further basic research indicates that “deficiency” is mainly reflected in the reduced number and diminished function of liver parenchymal cells and damage to the liver sinusoidal walls, while “blood stasis” is characterized by excessive ECM deposition and capillarization of liver sinusoids.9 In terms of treatment for liver-kidney yin deficiency, clinical manifestations include hypochondriac pain aggravated by exertion, soreness in the waist and knees, dry mouth and throat, restlessness, dizziness, insomnia, dry eyes, and a red tongue with scant white coating and minimal moisture. The pulse is string-like and fine. The therapeutic approach is to nourish the liver and kidney, with Yiguanjian (YGJ) commonly used in practice.10 YGJ, referenced in the Ancient Classical Famous Formulas Key Information Table (7 Formulas) and originating from the Yi Fang Jie Du, consists of ingredients such as Glehniae Radix (Bei Sha Shen), Ophiopogon tuber (Mai Dong), Angelica sinensis (Dang Gui), Rehmannia (Sheng Di Huang), Goji berry (Gou Qi Zi), and Melia toosendan (Chuan Lian Zi). However, as a complex Chinese herbal formula, YGJ’s components and targets are intricate, and the mechanism by which it treats liver fibrosis remains unclear.
Network pharmacology is a discipline that explores the relationships among drugs, targets, and diseases to virtually screen for potential active compounds, molecular targets, and signaling pathways, thus establishing “YGJ–target–disease” networks to predict the pharmacological basis and mechanisms of action of drugs.11 Recent technological advances have revolutionized this field, particularly with the advent of machine learning. The application of machine learning enables comprehensive analysis of extensive biomedical datasets, leading to the development of more accurate diagnostic methods and individualized patient treatments.12–14 Machine learning has enhanced our ability to discover disease biomarkers and elucidate underlying disease mechanisms.15 Mendelian randomization (MR), a method combining instrumental variables (IVs) and Mendelian inheritance laws, overcomes the limitations of traditional randomized controlled trials by minimizing confounding factors. With advancements in Mendelian theory and its growing practical application, MR is becoming an effective tool for inferring the causal relationship between drug targets (e.g., antagonists, agonists) and disease risk.16,17
In this study, we propose a comprehensive multilevel in silico validation strategy to elucidate the potential mechanisms of YGJ in treating liver fibrosis. This strategy leverages the complementary nature of multiple computational approaches: network pharmacology and machine learning were utilized for unbiased target screening and feature extraction; GEO expression analysis provided transcriptomic validation of disease relevance; MR offered genetic inference to establish putative causal relationships minimizing confounding biases; and molecular docking alongside dynamics simulations yielded biophysical evidence of direct YGJ–target interactions. By integrating these orthogonal layers of evidence, we aimed to systematically evaluate the active ingredients and potential targets of YGJ in the treatment of liver fibrosis. Furthermore, we performed gene expression and diagnostic validation, immune cell infiltration analysis, and single-gene GSEA enrichment analysis, with the detailed workflow depicted in Figure 1.
Materials and methods
Screening of core active ingredients and targets of TCM
The six herbal ingredients of YGJ (Sheng Di Huang, Bei Sha Shen, Mai Dong, Dang Gui, Gou Qi Zi, and Chuan Lian Zi) were input into the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP, https://old.tcmsp-e.com/tcmsp.php ) for screening. Active ingredients with oral bioavailability (OB) ≥ 30% and drug-like properties (DL) ≥ 0.18 were selected. The resulting data were merged and deduplicated to create a database of core active ingredients. These ingredients were then imported into Swiss Target Prediction (http://swisstargetprediction.ch/ ) for target prediction using a probability threshold of >0. The deduplicated target list was subsequently obtained. Additionally, active ingredients and targets were retrieved from the Encyclopedia of Traditional Chinese Medicine (ETCM) database (http://www.tcmip.cn/ETCM/index.php/Home/Index ). After further deduplication, the final YGJ-target database was compiled.
Identification of disease targets
Using “liver fibrosis” as the keyword, several relevant disease and control datasets were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/ ). These datasets were merged, and batch effects were removed using the Sangerbox online platform (http://sangerbox.com/ ), followed by differential analysis using the limma package (differentially expressed genes were identified using thresholds of |log2FC| ≥ 1.5 and adjusted P < 0.05) and weighted gene co-expression network analysis (WGCNA).18,19 The final disease target database was constructed by removing duplicate genes.
Immune cell infiltration analysis
Immune cell infiltration was estimated using the TIMER2.0 platform (http://timer.cistrome.org/ ) via the CIBERSORT algorithm, applied to the batch-corrected datasets from multiple sources.20 The CIBERSORT runtime parameters were set to their default values, and the sample filtering criterion was set to a CIBERSORT-calculated P-value < 0.05 to ensure the reliability of the estimated immune cell proportions.
Intersection of targets and construction of protein–protein interaction (PPI) network
The intersection of the active ingredient targets and disease targets was identified to obtain the core pharmacological targets for treating liver fibrosis. These intersection targets were imported into the STRING database (https://version-11-5.string-db.org/cgi/input?sessionId=bC4qb8Lym4Qw&input_page_show_search=on ), with “Homo sapiens” selected and a confidence threshold of >0.4. The original PPI network was generated, followed by optimization using Cytoscape 3.7.1. The MCODE plugin within Cytoscape 3.7.1 was employed to analyze the network, identifying subnetworks and feature gene sets based on the score values.21
Enrichment analysis
Enrichment analysis was conducted using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The intersection target list was imported into the Sangerbox online platform (http://sangerbox.com/ ) for analysis. GO enrichment analysis included three categories: biological process (BP), molecular function (MF), and cellular component (CC). Results were filtered with a significance threshold of P < 0.05.
Feature gene selection using machine learning algorithms
Machine learning techniques were employed to identify feature genes crucial to liver fibrosis. The dataset splitting strategy was as follows: the merged, batch-corrected dataset was randomly split into a training set and an internal validation set in a 7:3 ratio, and feature selection was performed only on the training set. To ensure robustness and comprehensive selection, Least Absolute Shrinkage and Selection Operator (LASSO) regression was used for variable selection and regularization, known for preventing overfitting by balancing model complexity and performance. The hyperparameter λ was selected using 10-fold cross-validation, and the optimal model was determined using the “1 standard error” criterion (λ = 1se) to identify genes with non-zero coefficients. The Random Forest (RF) algorithm, an ensemble learning method, was used to handle imbalanced data and estimate feature importance. It constructed 500 decision trees and evaluated feature importance using Recursive Feature Elimination (RFE), with the average decrease in the Gini coefficient serving as a measure of feature importance. Features with a score higher than 0.5 were considered significant. The XGBoost model, which integrates multiple weak classifiers (decision trees) to form a strong classifier, was also used. The hyperparameters for the XGBoost model were set as follows: the learning rate (eta) was 0.01, the maximum tree depth (max_depth) was 5, and the number of rounds (nrounds) was 500, with the optimal number of rounds determined via fivefold cross-validation. Feature importance was evaluated using the gain metric. Feature genes were identified from the intersection of the genes selected by LASSO logistic regression, XGBoost, and RF methods.
MR analysis
A two-sample MR analysis was conducted to investigate the causal relationship between feature genes and liver fibrosis. Single nucleotide polymorphisms (SNPs) were used as IVs. Cis-eQTL data for the target genes were obtained from the Integrated Epidemiology Unit (IEU) database (https://gwas.mrcieu.ac.uk/ ) and GTEx (http://www.gtexportal.org/ ) to serve as exposure factors. The specific data sources were as follows: cis-eQTL data for FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4 were derived from pooled whole-blood tissue data in the GTEx v8 database. Pooled GWAS data on liver fibrosis served as the outcome factor. The data were sourced from the 11th release (R11) of the Finnish database (FinnGen), with the dataset ID finn-b-K11_FIBROCHIRLIV. This dataset includes case and control samples of liver fibrosis and cirrhosis from individuals of European descent. SNPs with a P-value < 5 × 10−8 were selected, and IVs with strong linkage disequilibrium (LD) (r2 < 0.3, kb = 500) were excluded.22 MR analysis was conducted using the “TwoSampleMR” package (version 0.6.8) in R (version 4.4.1). Methods used included inverse variance weighted (IVW), MR-Egger, weighted median (WM), simple mode, and weighted mode. Sensitivity analyses were performed to assess the reliability of MR results using heterogeneity, horizontal pleiotropy, and leave-one-out (LOO) analysis. Heterogeneity was assessed, with P > 0.05 indicating no heterogeneity.
In silico gene expression cross-validation
To evaluate the generalizability of the selected signature genes, their differential expression and diagnostic performance were further validated on the training sets (GSE14323 and GSE164760) and independent validation sets (GSE6764 and GSE103580) to assess the model’s stability and generalizability across different data sources. Differential gene expression of biomarkers was validated in multiple independent datasets. Receiver operating characteristic (ROC) curves were generated to evaluate the diagnostic accuracy and specificity of five feature genes for liver fibrosis. Analysis and visualization were performed using the Sangerbox platform (http://sangerbox.com/login.html ).
Gene Set Enrichment Analysis (GSEA)
GSEA was performed using the Sangerbox platform (http://sangerbox.com/login.html ) and the KEGG database to identify significantly enriched pathways associated with the five feature genes and biomarkers in liver fibrosis. Samples were divided into high- and low-expression groups based on the median expression level of the feature gene, and 1,000 permutation tests were conducted, with a significance threshold set at |NES| ≥ 1.0 and FDR < 0.25.
Molecular docking
The 2D structures of active compounds were obtained from the PubChem database and stored in SDF format. The 3D structures of key targets were retrieved from the Protein Data Bank (PDB) in PDB format. Molecular docking was performed using Discovery Studio software, and the docking results were visualized with PyMOL. A redocking process was performed at the original active site to validate the docking reliability.
Molecular dynamics simulation
Molecular dynamics simulations were performed on the WeMol platform. The protein receptor was parameterized using the GROMACS 2023 software package and the CHARMM36 force field, and ligand topology files were generated using the CGenFF and GAFF force fields. The protein–ligand complex was placed in a TIP3P water box with 0.15 M NaCl to simulate physiological ionic strength. After energy minimization using the steepest descent method, a 100 ps NPT ensemble equilibrium was performed at 310 K (V-rescale thermostat) and 1 bar (Parrinello–Rahman pressure thermostat), followed by a 100 ns final molecular dynamics simulation. After the simulation, the trajectories were post-processed, and the root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) were calculated to assess the structural stability of the complex.
Results
Screening of core active ingredients and targets of TCM
A total of 59 active ingredients from five core herbs were identified from the TCMSP database after deduplication. These ingredients were then imported into the Swiss Target Prediction database, yielding 1,186 corresponding targets. After deduplication, 1,075 effective active ingredient–target pairs were obtained. Additionally, 735 targets were retrieved from the ETCM database (http://www.tcmip.cn/ETCM/index.php/Home/Index ). After removing duplicates, the final YGJ target database contained 1,198 unique targets.
Identification of disease targets
Four datasets (GSE14323, GSE164760, GSE6764, and GSE103580) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/ ), comprising a total of 260 samples from patients with chronic liver disease (covering various disease stages, including liver fibrosis, cirrhosis, and hepatocellular carcinoma) and 25 healthy control samples. Due to the lack of a unified disease staging standard in the original annotation information of each dataset and the limited sample size in some datasets, this study uniformly defined the above samples as the “disease group” to identify stable transcriptomic features associated with chronic liver disease status. It should be noted that the core objective of this study was to explore causal relationships between drug targets and diseases; subsequent MR analysis can provide causal inference evidence independent of sample composition. To mitigate batch effects, normalization was applied to the raw data of all datasets. The GSE14323 and GSE164760 datasets were combined to form the training cohort, while GSE6764 and GSE103580 were used as independent validation cohorts to confirm our findings. After batch effect removal and data merging via the Sangerbox platform (http://sangerbox.com/ ), the final dataset included 260 liver fibrosis samples and 25 control samples. Box plots revealed considerable distribution differences between the datasets before batch effect correction, indicating the presence of batch effects, whereas after batch effect removal, the data distributions aligned, with median values converging on the same line. Uniform Manifold Approximation and Projection (UMAP) analysis further confirmed the presence of batch effects before correction, with samples from each dataset clustering separately. However, after batch effect removal, samples from different datasets clustered together, indicating effective correction. Differentially expressed genes (DEGs) were identified using the limma package, with 144 upregulated genes and 54 downregulated genes. These results were visualized as volcano plots and heatmaps and organized into a matrix expression profile.
WGCNA was performed, and the soft threshold was determined to be 5 based on scale independence and average connectivity comparisons. Genes with similar expression profiles were classified into gene modules using hierarchical clustering based on TOM (topological overlap measure). The minimum module size for clustering was set to 30. Modules with a distance of less than 0.2 were merged, resulting in three co-expression modules. Notably, the black module, representing genes that could not be assigned to any module, was identified. Hub genes were filtered using the cutoff criterion of module membership (MM) > 0.3. After merging DEGs and hub genes and removing duplicates, the disease target database contained a total of 2,887 genes (see Fig. 2).
Immune cell infiltration
To explore the role of immune cells in liver fibrosis pathogenesis, immune cell infiltration analysis was performed to determine the proportions of 22 immune cell types in each sample. As shown in Figure 3a, correlation analysis revealed strong associations between certain immune cell types. For instance, naive B cells were strongly correlated with memory B cells (r = 0.33), effector B cells with monocytes (r = 0.32), and effector B cells with neutrophils (r = 0.41). Cytotoxic T lymphocytes exhibited strong correlations with resting CD4+ memory T cells (r = 0.58), follicular helper T cells (r = 0.42), and γδ T cells (r = 0.57). Resting CD4+ memory T cells were strongly correlated with follicular helper T cells (r = 0.53), regulatory T cells (r = 0.57), γδ T cells (r = 0.48), resting NK cells (r = 0.33), and M1 macrophages (r = 0.40). Further analysis revealed significant correlations between various immune cells, including regulatory T cells and γδ T cells (r = 0.48), regulatory T cells and resting NK cells (r = 0.43), and resting CD4+ memory T cells with neutrophils (r = 0.40), among others. Detailed correlations are shown in Figure 3b.
Box plot analysis further illustrated the differential expression of immune cells between the liver fibrosis and normal groups. Effector B cells, resting CD4+ memory T cells, γδ T cells, and M1 macrophages showed significantly higher expression levels in the liver fibrosis group than in the normal group (Fig. 3c).
Intersection of targets and construction of the PPI network
The intersection of the active ingredient–target library and the disease target library yielded 303 overlapping targets. These were imported into the STRING database to generate the original PPI network (Fig. 4). The PPI data were then analyzed using Cytoscape 3.7.1, calculating the degree, closeness, and betweenness values. The average degree value was 14.83, with closeness and betweenness values of 0.00127 and 507.64, respectively. Nodes with values lower than the average degree were excluded. Further analysis with the MCODE tool led to the identification of a subnetwork (Fig. 4).
Enrichment analysis
A total of 2,420 entries were obtained from GO functional analysis. These included 1,965 BP entries, 193 CC entries, and 262 MF entries. The top 10 enriched terms, based on the number of enriched genes, were visualized (Fig. 5). KEGG pathway analysis identified 148 pathways, from which the top 20 were selected based on both the number of enriched genes and P-values. These pathways were primarily associated with non-alcoholic fatty liver disease, AMPK, AGE–RAGE signaling, and tyrosine metabolism, among others. The enrichment bubble chart and bar chart illustrating these pathways are shown in Figure 5.
Feature gene selection using machine learning algorithms
To identify the 303 key genes, we applied various machine learning algorithms. Specifically, LASSO regression identified 23 genes, yielding a ROC curve with an area under the curve (AUC) of 0.96. The XGBoost algorithm identified 10 genes, while the RF method identified 47 important genes. A Venn diagram was used to identify common genes across all three methods, resulting in the selection of five feature genes, as shown in Figure 6.
MR analysis
MR analysis was performed to confirm the causal relationship between key target genes and liver fibrosis using methods such as IVW, MR-Egger, WM, simple mode, and weighted mode. The results, displayed in a forest plot, indicated that five feature genes (FABP4 (GTEx Database), MDM2 (eQTL-a-ENSG00000135679), AKR1B1 (eQTL-a-ENSG00000085662), PDGFRB (eQTL-a-ENSG00000085662), and NR1H4 (GTEx Database)) exhibited strong causal associations with liver fibrosis (odds ratio (OR) > 0.9), suggesting their potential as effective drug targets. Among these, the ORs for FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4 exceeded 1, indicating that their expression levels are positively correlated with the likelihood of liver fibrosis and classifying them as risk factors. Consequently, we hypothesize that the active compounds in YGJ may exert their therapeutic effects by acting as antagonists or downregulators of these risk-associated genes. Sensitivity analyses, including funnel plot evaluation, confirmed the robustness of these MR results. Leave-one-out analysis demonstrated consistent MR findings after systematically removing each SNP, further supporting the reliability of these results (Fig. 7).
Gene expression validation
The expression of the five key genes (FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4) was evaluated in both the training (GSE14323 and GSE164760) and independent validation cohorts. The expression trends for these genes were consistent across all datasets (Fig. 8). Specifically, compared with the control group, the expression levels of FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4 were higher in the liver fibrosis group, highlighting these core targets as potential risk factors. Furthermore, ROC curve analysis revealed AUC values greater than 0.5 for all five genes, confirming their potential diagnostic ability for liver fibrosis (Fig. 9). These results indicate that FABP4, AKR1B1, PDGFRB, and NR1H4 exhibit consistent differential expression in chronic liver diseases, suggesting their potential as biomarkers for chronic liver diseases. It is worth noting that, since the case cohort includes samples from different disease stages, the aforementioned diagnostic performance requires further validation in a more strictly defined cohort of early-stage liver fibrosis. It is also worth noting that MDM2 did not demonstrate outstanding diagnostic performance (with an AUC below 0.5) in the training set or in GSE82107. This may be due to differences in the tissues in which the same gene is expressed across GEO datasets; its diagnostic performance requires further research and validation.
GSEA
To elucidate the underlying mechanisms, GSEA was conducted to analyze the significantly enriched pathways for the five core targets (FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4). FABP4 was mainly enriched in actin cytoskeleton regulation pathways; MDM2 in apoptosis and Wnt signaling pathways; AKR1B1 in PPAR and pyruvate metabolism pathways; PDGFRB in ECM–receptor interaction and MAPK signaling pathways; and NR1H4 in alanine, aspartate, and glutamate metabolism and tryptophan metabolism pathways (Fig. 10).
Molecular docking and binding mode analysis
The top 10 effective ingredients from the “Core TCM–Active Ingredients–Disease Pathways–Intersecting Targets” network, based on docking scores, were selected as core compounds: Cytidylic Acid B, Adenosine, Adenine Nucleoside, Stigmasterol, Cytidylic Acid A, Quercetin, Uridylic Acid, Bergaptin, Cnidilin, and Ethyl Linolenate. Molecular docking was performed with FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4. The highest binding scores were obtained for the following core compound–target pairs: AKR1B1–Stigmasterol, AKR1B1–Bergaptin, AKR1B1–Melianone, and NR1H4–Melianone, with average binding energies of −14.64, −12.86, −16.62, and −12.52 kcal/mol, respectively. Binding energies below −4.25, −5.0, or −7.0 kcal/mol indicate moderate, good, or strong binding affinity, respectively, suggesting that all four core compound–target pairs exhibit strong binding activity (Fig. 11).
Molecular dynamics simulation
To verify binding stability, molecular dynamics (MD) simulations were performed on the docking complexes over a 100 ns period with a time step of 0.002 ns. RMSD was used to evaluate the stability of protein–ligand conformations, while RMSF was employed to analyze positional variations of proteins under specific pressure and temperature conditions. The results indicated that the RMSD curves for all pairs displayed good stability over the 100 ns simulation period. Specifically, AKR1B1–Stigmasterol stabilized after 25 ns, AKR1B1–Bergaptin after 20 ns, AKR1B1–Melianone after 10 ns, and NR1H4–Melianone after 20 ns. RMSF peaks indicated significant fluctuations in specific regions throughout the simulation. The RMSF values for the four core compound–target pairs ranged from 0.01 Å to 0.69 Å, suggesting that these protein–ligand complexes remained stable under the simulated conditions. Overall, these molecular dynamics simulations support the stability of the selected docking complexes (Fig. 12).
Discussion
Through PPI network construction and machine learning-based screening, five core targets (FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4) were computationally predicted as putative risk factors for liver fibrosis. Given their identification as risk factors in our MR analysis (OR > 1), we hypothesize that YGJ potentially exerts its antifibrotic effects by inhibiting or downregulating these specific targets. It should be noted that these findings are derived from bioinformatic inferences and require further experimental corroboration. FABP4 is a lipid factor secreted by adipocytes and macrophages, and studies have demonstrated its involvement in lipid metabolism and inflammatory responses.23 MDM2, an E3 ubiquitin ligase, plays a crucial role in regulating the degradation of various proteins, including p53/TP53, which is implicated in hepatic stellate cell-mediated ferroptosis during liver fibrosis.24,25 AKR1B1 encodes an aldose reductase enzyme, which is involved in cancer metabolism and intracellular fructose supply.26 PDGFRB, a receptor tyrosine kinase, activates key signaling pathways such as PI3K/Akt, RAS/MAPK, and JAK/STAT, promoting cell proliferation and survival.27 NR1H4, which encodes a bile acid receptor, is crucial for liver development and has been linked to cholestatic liver diseases.28
The findings suggest that these five key genes play distinct roles in the pathological progression of liver fibrosis. Furthermore, the “active ingredient–target network” constructed in this study indicates that YGJ exerts its therapeutic effects not through a single target or active ingredient but via a “multi-ingredient, multi-target” mechanism. This highlights the complex interactions between different active ingredients and targets in the treatment of liver fibrosis.
Immune infiltration analysis revealed an increased presence of effector B cells, resting CD4+ memory T cells, γδ T cells, and M1 macrophages during the development of liver fibrosis, suggesting that immune cells play a pivotal role in disease progression.
KEGG and GSEA enrichment analyses demonstrated that the key targets of YGJ in the treatment of liver fibrosis are enriched in multiple signaling pathways, including MAPK, AGE–RAGE, tyrosine metabolism, actin cytoskeleton regulation, apoptosis, Wnt signaling, PPAR signaling, pyruvate metabolism, ECM–receptor interactions, alanine, aspartate and glutamate metabolism, and tryptophan metabolism. MAPK signaling can induce mitochondrial-dependent apoptosis, promoting Bax translocation to mitochondria, inhibiting Bcl-2, and activating caspase-3, which may work synergistically with the Bcl-2/Bax signaling pathway in regulating liver fibrosis.29 The AGE–RAGE signaling pathway has been identified as a key mechanism in the activation of hepatic stellate cells.30 Wnt signaling regulates lactate dehydrogenase-A nuclear translocation, mediating glycolysis in hepatic stellate cells during the development of liver fibrosis.31 Targeting PPAR-γ has been shown to inhibit hepatic stellate cell activation and improve liver fibrosis.32 The role of pyruvate metabolism in driving hepatic stellate cell activation depends on its oxidative phosphorylation in the mitochondria or its metabolic processing in the cytoplasm.33
Furthermore, MR was employed to explore the potential causal relationships between the treatment targets (FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4) and liver fibrosis. This approach leverages eQTL and liver fibrosis association data to establish causal links, thereby enhancing the predictive validity of the identified targets. MR analysis, employing five different methods to estimate causal effects, combined with sensitivity analysis and colocalization verification, provides robust causal estimates while minimizing the potential influence of pleiotropy in the IVs, further supporting the relevance of these targets in liver fibrosis.
In MR analysis, the five feature genes (FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4) exhibited strong causal relationships with liver fibrosis (OR > 0.9), and all demonstrated significant potential as drug targets. OR values greater than 1 suggest that increased expression of these genes correlates with a higher likelihood of liver fibrosis, positioning them as risk factors for the disease. From a pharmacological perspective, this directionality suggests that therapeutic interventions, such as YGJ, would need to inhibit or suppress these targets to mitigate disease progression. This aligns with the observed expression differences of these genes in liver fibrosis patients compared with control groups across all datasets. The consistency between MR analysis and differential expression analysis further corroborates the robust causal relationship of FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4 with liver fibrosis. ROC curves across all datasets confirm the diagnostic potential of these genes, with AUC values consistently exceeding 0.5, demonstrating predictive performance. At the same time, a major limitation of this study is that the case cohort in the GEO dataset encompasses a continuous disease spectrum ranging from liver fibrosis to hepatocellular carcinoma, and precise stratification based on existing annotation information is not possible. However, the core findings of this study—namely, the causal relationships between FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4 and liver disease risk—are primarily based on MR analysis. MR analysis utilizes genetic variants as IVs, and its causal estimates are not influenced by the composition of the case cohort (since the association between exposure and outcome is based on population-level genetic data rather than transcriptomic data from the case cohort). Therefore, despite phenotypic heterogeneity in the transcriptomic-level validation, the MR results provide robust support for the key causal inferences of this study.
Lastly, molecular docking and molecular dynamics simulations were employed to validate the binding affinity of key active ingredients in YGJ with the identified core targets. The binding energies between the core ingredients and core targets suggest strong binding interactions, providing further evidence of the synergistic therapeutic effects of these ingredients in treating liver fibrosis. However, further research is needed to explore whether these interactions extend to other active ingredients in YGJ.
In summary, our study successfully integrates network pharmacology and MR to establish a causal relationship between “YGJ–key components–core targets–liver fibrosis,” identifying and validating FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4 as reliable biomarkers for liver fibrosis with good predictive performance. Functional investigations reveal significant interactions between these targets and immune-related pathways in liver fibrosis. These findings offer valuable insights into the complex molecular mechanisms of liver fibrosis and emphasize the potential for developing personalized diagnostic and therapeutic strategies, paving the way for targeted and effective interventions.
This study has several limitations:
Due to the characteristics of TCM in syndrome differentiation and treatment, the dosage ratio of YGJ for treating various types of liver fibrosis (primarily targeting Liver-Kidney Yin Deficiency Syndrome) may vary, potentially influencing its efficacy. Insufficient clinical data were collected to further refine the liver fibrosis model and perform weighted analyses of different syndromes and dosages. Furthermore, Chinese herbal formulations involve complex “multi-ingredient, multi-target” actions, and quantifying synergistic interactions between active ingredients and their activation or inhibition of targets requires advanced network algorithms, which remains an area for future improvement in network pharmacology. With the advancement of artificial intelligence, network pharmacology, and molecular docking technologies, TCM and its formulations are likely to be more effectively integrated with tissue engineering for clinical application.
In the MR analysis of drug targets, stringent criteria were used; therefore, other SNPs with causal relationships to liver fibrosis may still be unidentified. Due to sample limitations, stratification by disease type and patient sex was not performed, which may have introduced some bias. Further clinical validation of these key targets will require additional datasets.34
The GEO dataset used in this study may present variability in gene expression across different tissues, which could potentially affect AUC values and other results.
The GEO dataset used for transcriptomic validation includes samples of chronic liver disease at various stages of progression; however, due to limitations in annotation of the raw data, it was not possible to perform precise subgroup analyses or exclude samples with hepatocellular carcinoma. Consequently, ROC analysis results may be influenced by disease severity and cannot be directly extrapolated to the diagnostic context of early-stage liver fibrosis. However, the main conclusion of this study—the causal relationship between core genes and liver disease risk—is supported by MR analysis, which is independent of the aforementioned sample composition. Independent validation in a strictly defined early-stage fibrosis cohort remains necessary in the future.
Finally, it is important to emphasize that this is primarily a computational and genetic inference-based study. Although multiple algorithms and independent datasets were used for cross-validation, the precise mechanistic pathways lack direct in vivo or in vitro experimental validation. Therefore, the findings should be interpreted as predictive hypotheses intended to guide future experimental research. Comprehensive in vitro and in vivo studies are still required to empirically confirm these pharmacological targets and pathways.
Conclusions
Network pharmacology and MR identified FABP4, MDM2, AKR1B1, PDGFRB, and NR1H4 as potential risk biomarkers for liver fibrosis, suggesting that YGJ may exert therapeutic effects by downregulating or antagonizing these targets. Functional analyses further revealed significant interactions between these targets and immune-related pathways. These findings offer mechanistic insights and potential strategies for targeted intervention, although direct experimental validation is still needed to confirm these in silico predictions.
Declarations
Acknowledgement
We thank the Guangxi University of Chinese Medicine (Guangxi, China) for providing laboratory equipment and technical support.
Ethical statement
Not applicable.
Data sharing statement
The data used and analyzed in this study are included within the article.
Funding
This study was supported by the National Natural Science Foundation of China (Grant No. 82204755), the Key Laboratory of Faculty of Chinese Medicine Science, Guangxi University of Chinese Medicine (No. 02: Basic Research on the Prevention and Treatment of Chronic Liver Diseases with Chinese Medicine and Zhuang-Yao Medicine), the Guangxi Natural Science Foundation (Grant Nos. 2023GXNSFBA026274 and 2024GXNSFAA010235), the Guangxi Zhuangyao Pharmaceutical Key Laboratory (Grant No. GXZYYKF2023-05), the Guangxi University of Traditional Chinese Medicine School-level Project Youth Fund (Grant No. 2022QN008), the Faculty of Chinese Medicine Science, Guangxi University of Chinese Medicine Research Project (Grant Nos. 2024ZZA002, 2024ZZA003, 2024ZZB006, 2024ZZB007, and 2024ZZB010), and the National Student Innovation and Entrepreneurship Training Programme of the Faculty of Chinese Medicine Science at the Guangxi University of Chinese Medicine (Grant No. 202513643017).
Conflict of interest
The authors declare that there are no conflicts of interest.
Authors’ contributions
Study concept and design (YKJ, JHW), acquisition of data (LW), analysis and interpretation of data (JHW, YZ), drafting of the manuscript (YKJ), critical revision of the manuscript for important intellectual content (TJZ), administrative, technical, or material support (HYX), and study supervision (RWZ, HYX). All authors have made significant contributions to this study and have approved the final manuscript.