Integrated transcriptomic and co-expression network analysis identifies HIF1A as a key immune regulator in biliary atresia

Introduction

Biliary atresia (BA) is a progressive cholangiopathy unique to the neonatal period, characterized by an initial insult to the extrahepatic bile ducts, which subsequently leads to bile duct obstruction, fibrosis, and ultimately results in severe cholestasis, hepatic cirrhosis, and end-stage liver failure (1). Kasai portoenterostomy (KPE) remains the first-line surgical intervention and can partially restore bile flow if performed early. However, data from multicenter cohorts reveal that the native liver survival (NLS) rate within 5 years post-KPE is only 30–50%, with merely 15–20% of patients maintaining native liver function into adulthood (2). The majority of BA patients ultimately require liver transplantation during childhood or adolescence, underscoring the urgent need to elucidate the underlying pathogenesis of BA and to improve strategies for early diagnosis and targeted therapeutic intervention (3,4).

The etiology and pathogenesis of BA are multifactorial and remain incompletely understood. Perinatal viral infections, such as rotavirus and cytomegalovirus (CMV)—alongside congenital genetic susceptibility, immune dysregulation, and environmental factors have all been implicated as potential contributors to disease initiation and progression (5-7). Histopathological analyses of liver tissue from BA patients reveal extensive infiltration by T lymphocytes, macrophages, and neutrophils, accompanied by marked upregulation of chemokines and proinflammatory cytokines. This persistent proinflammatory immune microenvironment exacerbates cholangiocyte injury and drives the progression of hepatic fibrosis (8,9).

In addition, signaling pathways critical for embryonic development play pivotal roles in BA pathogenesis. Aberrant activation of the Hedgehog signaling pathway promotes excessive hepatic progenitor cell proliferation and ductular reactions, and facilitates fibrogenesis through transforming growth factor-β (TGF-β)-mediated pathways (10,11). Dysregulation of the NOTCH signaling pathway further impairs proper bile duct cell differentiation and contributes to fibrotic remodeling (12). Moreover, hypoxia-inducible factor-1α (HIF-1α) is significantly upregulated in biliary epithelial cells of BA patients and may promote fibrotic and inflammatory gene expression through cross-talk with the NF-κB signaling axis, suggesting its potential as a novel therapeutic target (13).

In recent years, the rapid advancement of high-throughput transcriptomic sequencing and systematic bioinformatics approaches has enabled researchers to investigate gene expression alterations, dysregulated signaling pathways, and critical gene modules associated with the pathogenesis of BA from a systems-level perspective. By applying strategies such as differential gene expression analysis, functional enrichment, protein-protein interaction (PPI) network construction, and weighted gene co-expression network analysis (WGCNA), multiple candidate genes and regulatory modules have been identified that are associated with BA molecular subtypes and clinical outcomes (14-16). These findings offer novel insights into the development of early diagnostic biomarkers and precision therapeutic targets.

In this study, we systematically integrated multiple BA-related transcriptomic datasets from the Gene Expression Omnibus (GEO) database to identify and validate key developmental- and immune-related genes. We further constructed a core gene interaction network to elucidate the molecular mechanisms underlying BA pathogenesis, thereby providing a theoretical basis for future experimental validation and clinical intervention strategies. We present this article in accordance with the STREGA reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-590/rc).

Methods

GEO data acquisition

BA-related microarray gene expression datasets GSE46960, GSE159720, and GSE15235 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi). Detailed information regarding each dataset is summarized in Table 1. In GSE46960, 14 control liver tissues were included as non-BA pediatric samples obtained during unrelated surgical procedures and confirmed to be free of cholestasis or inflammation according to the original dataset annotation. Detailed demographic information, such as age matching, was not uniformly available and thus could not be fully adjusted.

Among the three datasets, GSE46960 served as the training set for differential expression and co-expression analyses. GSE159720 and GSE15235 were used for external validation and extended functional analysis. Despite its small sample size (three controls, four BA cases), GSE159720 was generated on a different microarray platform, offering independent validation of expression trends. GSE15235, though lacking normal controls, contains many BA samples of two subtypes (fibrotic and inflammatory) and was used to assess intra-disease variation in HIF1A expression and conduct downstream functional and immune analyses based on HIF1A stratification. Comprehensive differential expression statistics for the training cohort are summarized in table available at https://cdn.amegroups.cn/static/public/tp-2025-590-1.xlsx, and effect sizes with 95% confidence intervals (CIs) for the validation datasets are provided in tables available at https://cdn.amegroups.cn/static/public/tp-2025-590-2.xlsx.

All included datasets were derived from previously published studies that had obtained local institutional review board approval and informed consent according to the original reports. This secondary analysis used only de-identified, publicly available data and did not involve any direct patient contact or additional sample collection. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Data preprocessing

The raw microarray files (e.g., CEL for Affymetrix; IDAT or other vendor-native formats as applicable) for GSE46960 (GPL6244), GSE159720 (GPL18573), and GSE15235 (GPL570) are available in GEO. This analysis used the submitter-provided normalized expression matrices; platform-appropriate preprocessing is documented in the corresponding GEO records. Probe IDs were mapped to official gene symbols using platform-specific Bioconductor annotation packages, and when multiple probe sets mapped to the same gene, expression was summarized by the median (ties were retained where identical values occurred).

Array-level quality control (QC) metrics [e.g., normalized unique solution enrichment (NUSE)/relative log error (RLE)] and array weights were not recomputed from raw intensities; accordingly, no array-level exclusion was performed based on NUSE/RLE thresholds. We present principal component analysis (PCA) of normalized matrices as an exploratory view, and platform-wise analyses were kept separate to avoid cross-platform artifacts. Recomputing QC from heterogeneous raw files might introduce additional batch effects and was therefore not performed.

GSE46960 served as the training dataset, whereas GSE159720 and GSE15235 were exploratory validation cohorts. GSE159720 (three controls, four BA) provided directional confirmation, while GSE15235, though control-free, enabled within-disease analysis between fibrotic and inflammatory subtypes.

Clinical variables (e.g., age, CMV status, fibrosis stage) were unavailable; thus, ComBat or SVA correction was not applied to prevent over-adjustment. Analyses emphasized directional consistency and effect-size interpretation [log₂fold change (log₂FC), 95% CI] over nominal P values, consistent with best practice for small or control-free transcriptomic datasets (17,18). A PCA overview is provided in Figure S1; no per-array QC table is included. Effect sizes with 95% CIs have been summarized in tables available at https://cdn.amegroups.cn/static/public/tp-2025-590-2.xlsx as described above.

Differential gene expression analysis

Differentially expressed genes (DEGs) were identified with the limma package (R), modeling BA vs. control with the control group as the reference. To control for multiple testing at the gene level, we applied the Benjamini-Hochberg false discovery rate (BH-FDR) to limma’s moderated statistics. Genes with q<0.05 (BH-FDR) and |log₂FC| >0.58 (~1.5-fold) were considered significant; the corresponding q-values are reported in table available at https://cdn.amegroups.cn/static/public/tp-2025-590-1.xlsx. This FC threshold balances sensitivity and biological relevance and is commonly used in transcriptomic analyses (19,20). Full DEG statistics are provided in table available at https://cdn.amegroups.cn/static/public/tp-2025-590-1.xlsx.

To assess robustness, we repeated the DE analysis under alternative |log₂FC| cutoffs (0.3 and 1.0) with FDR <0.05. The direction of effects was preserved, and key signals (including HIF1A and immune-related genes) persisted under these alternatives.

For downstream exploratory procedures [Gene Ontology (GO)/Kyoto Encyclopedia of Genes and Genomes (KEGG)/Gene Set Enrichment Analysis (GSEA) enrichment, gene correlations, and WGCNA], we generally retained nominal P<0.05 given the limited number of planned comparisons and small sample sizes. For families of immune cell-type comparisons, BH-FDR was applied where applicable. Limiting FDR control to large-scale gene-level testing while reporting nominal P for targeted exploratory analyses follows best practice for hypothesis-generating studies in small public cohorts (21,22).

For visualization, a volcano plot and a top-gene heatmap were produced; code to reproduce these figures (and any previews made on hosted tools) is provided in the Supplementary scripts. To explore putative functions and pathways, GO and KEGG enrichment analyses were conducted with clusterProfiler under the above analysis policy.

Screening and validation of BA development-related differential genes

Relevant gene sets associated with the Hedgehog and NOTCH signaling pathways in BA were retrieved from the GSEA database (https://www.gsea-msigdb.org/gsea/index.jsp) following a comprehensive literature search on PubMed. DEGs identified from the training set (GSE46960) were intersected with these pathway gene sets, and the overlap was visualized using a Venn diagram. The expression profiles of the intersected genes were subsequently validated in the independent validation set (GSE159720).

The Hedgehog (HALLMARK_HEDGEHOG_SIGNALING, M5890) and NOTCH (HALLMARK_NOTCH_SIGNALING, M5891) gene sets were retrieved from MSigDB v7.5.1 Hallmark (M5890/M5891) (23). HIF1A and RPS27A are recognized members of the NOTCH-related regulatory network according to MSigDB and Reactome annotations, where HIF1A participates in NOTCH-mediated transcriptional activation under hypoxia, and RPS27A is involved in ubiquitin-dependent signal modulation. Because this study aimed to prioritize shared development-related signals rather than quantify pathway activity, we applied binary intersection to identify common genes between DEGs and pathway members. Nevertheless, we acknowledge that gene set variation analysis (GSVA) or single-sample GSEA (ssGSEA) could provide continuous pathway activity scores, and this approach will be considered in future analyses.

Screening of immune-related differential genes in BA

Immune-related DEGs were obtained by intersecting DEGs from the training set (GSE46960) with gene sets associated with 28 types of immune infiltrating cells. A heatmap of these immune-related DEGs was generated using the “heatmap” package in R. PPI analysis was conducted by submitting the relevant genes to the STRING database (https://cn.string-db.org/). The resulting PPI network was visualized and further analyzed using Cytoscape software (version 3.8.0). Additionally, a correlation heatmap of the selected genes was generated using the Bioinformatics online platform (https://www.bioinformatics.com.cn/).

Differential gene screening based on HIF1A expression

Patients in the GSE46960 were stratified into high and low groups based on the median expression level of HIF1A. Differential gene expression analysis between the two groups was performed using the “limma” package in R. Volcano plots and heatmaps of the top 30 DEGs were generated via the Hiplot online platform. Functional enrichment analyses, including GO and KEGG pathways, were conducted using the “clusterProfiler” package. Additionally, GSEA was performed with GSEA software version 4.3.0, utilizing the c2 gene set from the Molecular Signatures Database (MSigDB).

Subtyping analysis of the GSE15235 data set

Patients in the validation set (GSE15235) were classified into inflammatory and fibrotic subtypes, and differences in HIF1A expression between these groups were analyzed. Subsequently, inflammatory subtype patients were stratified based on HIF1A expression levels for differential gene expression analysis. DEGs from the training set (GSE46960) and validation set (GSE15235) were intersected with immune-related gene sets, and the overlaps were visualized using Venn diagrams. PPI analyses and correlation assessments were performed, with interaction networks constructed using Cytoscape. Correlation coefficients were calculated via Sangerbox software, and corresponding correlation heatmaps were generated using the Bioinformatics online platform.

Immune infiltration analysis

Immune-cell enrichment was performed using xCell (default settings) (24). Because xCell and LM22/CIBERSORT enumerate non-identical and only partially overlapping cell categories, cross-method validation was not undertaken to avoid introducing additional assumptions. xCell results are reported as exploratory.

Correlation analysis between HIF1A expression and key immune indicators, including neutrophil/macrophage markers and major histocompatibility complex (MHC) molecule expression, was performed using Sangerbox and Bioinformatics.com.cn platforms, and visualized as correlation heatmaps.

BH-FDR correction was applied where applicable (e.g., multiple cell-type comparisons) to control FDRs, while group-level comparisons retained nominal P<0.05 due to small sample sizes and limited independent tests. This combined approach balances statistical rigor and interpretability for small-sample transcriptomic analyses.

WGCNA

WGCNA was conducted using the WGCNA package in R to explore gene co-expression patterns associated with HIF1A. DEGs identified between high- and low-HIF1A expression groups within the inflammatory subtype of the validation dataset (GSE15235) were used as input. This DEG-based WGCNA approach was designed to focus on genes transcriptionally linked to HIF1A and enhance the interpretability of co-expression relationships specifically relevant to HIF1A-driven regulatory programs.

The soft-threshold power (β) was determined as seven using the pickSoftThreshold function, achieving a scale-free topology fit index (R²) of 0.89 (Figure S2). The minimum module size was set to 10, mergeCutHeight to 0.25, and deepSplit to 2 to balance network granularity and stability. Modules with dissimilarity lower than 0.25 were merged automatically. Module stability was verified by 1,000 random permutations and eigengene preservation analysis, confirming robust clustering consistency {diagnostics in Figure S3A-S3C: scale-free topology fit vs. β, mean connectivity vs. β, and sensitivity across β ∈ [6, 9]}.

Modules exhibiting the highest correlation with HIF1A expression were selected for further analysis. Key genes within the most correlated modules were extracted and subjected to functional enrichment analysis using the Metascape platform to characterize biological functions and pathways.

Statistical analysis

All statistical analyses were performed using R software (version 4.3.0; R Foundation for Statistical Computing, Vienna, Austria) and GraphPad Prism (version 9.0; GraphPad Software, San Diego, CA, USA), unless otherwise specified. For two-group comparisons of continuous variables (e.g., gene expression levels, xCell enrichment scores), Student’s t-test was applied for approximately normally distributed data; otherwise, non-parametric tests (e.g., Wilcoxon rank-sum test) were considered in sensitivity analyses. For multiple-group comparisons, one-way analysis of variance (ANOVA) with appropriate post-hoc tests or Kruskal-Wallis tests was used where applicable.

Differential expression analysis was conducted with the limma package using empirical Bayes moderation, and P values were adjusted for multiple testing using the FDR method at the gene level, with q<0.05 and |log₂FC| >0.58 as primary thresholds. For downstream enrichment, correlation, and network analyses, nominal P<0.05 was generally considered statistically significant, and FDR correction was applied selectively where families of comparisons were large, as detailed in the corresponding subsections. All tests were two-sided.

Results

Differential expression analysis and GO/KEGG enrichment

Using the Control group as the reference, a total of 343 DEGs were identified (BH-FDR q<0.05 and |log₂FC| >0.58) between BA and Control groups, including 295 upregulated and 48 downregulated genes (Figure 1A). The heatmap in Figure 1B illustrates the expression profiles of the top 15 significantly upregulated and downregulated DEGs. Notably, genes such as SPP1, EPCAM, and TIMP1 exhibited elevated expression levels in BA samples, whereas SLC4A1, KCNMA1, and ITGA2B showed decreased expression. Comparable DEG patterns were observed when alternative log₂FC thresholds (0.3 and 1.0) were applied, confirming the robustness of expression trends (see table available at https://cdn.amegroups.cn/static/public/tp-2025-590-1.xlsx for the complete gene list).

Figure 1 Identification and functional enrichment analysis of DEGs in BA. (A) Volcano plot displaying a total of 343 DEGs between the BA and control groups, including 295 upregulated genes (red) and 48 downregulated genes (blue). (B) Heatmap showing the top 15 most significantly upregulated and downregulated DEGs. (C) GO enrichment analysis of upregulated genes. (D) KEGG pathway enrichment analysis of upregulated genes, with significant enrichment in the TNF signaling pathway, cytokine-cytokine receptor interaction, ECM-receptor interaction, and focal adhesion pathways. (E) GO enrichment analysis of downregulated genes. (F) KEGG enrichment analysis of downregulated genes, showing significant enrichment only in the malaria-related pathway. BA, biliary atresia; BP, biological process; CC, cellular component; DEG, differentially expressed gene; ECM, extracellular matrix; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; N, normal; T, tumor; TNF, tumor necrosis factor.

GO enrichment analysis of upregulated genes revealed significant enrichment in biological processes (BPs) related to chemotaxis, chemokine-mediated signaling, and wound healing; cellular components (CCs) including cell-matrix junction, focal adhesion, and collagen-containing extracellular matrix (ECM); as well as molecular functions (MFs) such as glycosaminoglycan binding and cytokine activity (Figure 1C). KEGG pathway analysis demonstrated that upregulated genes were significantly enriched in pathways including the tumor necrosis factor (TNF) signaling pathway, cytokine-cytokine receptor interaction, ECM-receptor interaction, focal adhesion, and viral protein interactions with cytokines (Figure 1D).

For downregulated genes, GO analysis indicated enrichment in BP associated with erythrocyte development, one-carbon compound transport, and gas transport; CC terms such as hemoglobin complex and basolateral plasma membrane; and MF terms including cytoskeletal structural constituent and oxygen binding (Figure 1E). KEGG pathway analysis of downregulated genes showed significant enrichment only in the malaria-related pathway (Figure 1F).

Identification of development and immune-related DEGs

Literature mining through the PubMed database identified the Hedgehog and NOTCH signaling pathways as closely associated with the pathogenesis of BA. Gene sets for these pathways were obtained from the GSEA database, including 205 genes for the Hedgehog pathway and 253 genes for the NOTCH pathway. Intersecting the DEGs from the training data set (GSE46960) with the Hedgehog pathway genes yielded two development-related DEGs: TGFB2 and RPS27A (Figure 2A). Intersection with the NOTCH pathway genes identified five development-related DEGs: HIF1A, ELF3, ACTA2, RPS27A, and MYC (Figure 2B). After deduplication, a total of six development-related genes were further validated in the validation data set (GSE159720), where only HIF1A showed a statistically significant difference in expression between BA and control groups (Figure 2C).

Figure 2 Identification and validation of development and immune-related DEGs in BA. (A) Venn diagram showing the overlap between DEGs from GSE46960 and Hedgehog pathway genes, identifying two development-related DEGs: TGFB2 and RPS27A. (B) Venn diagram showing the overlap between DEGs from GSE46960 and NOTCH pathway genes, identifying five development-related DEGs: HIF1A, ELF3, ACTA2, RPS27A, and MYC. (C) Box plots showing expression of six development-related DEGs in the validation dataset GSE159720; only HIF1A showed a statistically significant difference between BA and control groups. (D) Venn diagram showing the intersection of DEGs from GSE46960 and immune cell infiltration-related genes, identifying 29 immune-related DEGs. (E) Heatmap showing expression profiles of HIF1A and the 29 immune-related DEGs in the GSE46960 dataset. (F) PPI network constructed using the STRING database, indicating HIF1A as a central hub gene among the 29 immune-related DEGs. (G) Correlation analysis demonstrating that HIF1A expression is positively correlated with the majority of immune-related DEGs. ***, P<0.001. BA, biliary atresia; DEG, differentially expressed gene; N, normal; NC, normal control; PPI, protein-protein interaction; T, tumor.

In addition, intersection of the 343 DEGs from GSE46960 with 782 immune cell infiltration-related genes resulted in 29 immune-related DEGs (Figure 2D). Heatmap analysis revealed distinct expression differences between BA and control groups for these immune-related genes (Figure 2E). A PPI network was then constructed using the STRING database, where HIF1A was located at the center, potentially interacting with multiple immune-related genes (Figure 2F). Correlation analysis further demonstrated that HIF1A expression was positively correlated with most of these immune-related genes (Figure 2G).

Differential gene expression and functional enrichment analysis based on HIF1A expression

Patients in the training dataset GSE46960 were stratified into high and low HIF1A expression groups based on the median expression level of HIF1A. A total of 461 DEGs were identified, including 382 upregulated and 79 downregulated genes (Figure 3A). The top 15 upregulated and downregulated DEGs are illustrated in the heatmap (Figure 3B).

Figure 3 Differential gene expression, functional enrichment, and GSEA analysis based on HIF1A expression levels. (A) Volcano plot showing a total of 461 DEGs between high and low HIF1A expression groups, including 382 upregulated genes (red) and 79 downregulated genes (blue). (B) Heatmap displaying the expression profiles of the top 15 upregulated and downregulated DEGs between the high and low HIF1A groups. (C) GO enrichment analysis of upregulated genes reveals significant involvement in response to bacterial-derived molecules, response to lipopolysaccharide, regulation of inflammatory responses, and cytokine production. (D) KEGG pathway analysis of upregulated genes shows enrichment in TNF signaling, cytokine-cytokine receptor interaction, HIF-1 signaling, and other inflammation- and immunity-related pathways. (E) GO enrichment analysis of downregulated genes indicates roles in fatty acid metabolism, biosynthetic processes, and transport functions. (F) KEGG analysis of downregulated genes demonstrates enrichment in ABC transporters, PPAR signaling, insulin resistance, and bile acid metabolism pathways. ABC, ATP-binding cassette; ATP, adenosine triphosphate; BP, biological process; CC, cellular component; DEG, differentially expressed gene; FC, fold change; GO, Gene Ontology; GSEA, Gene Set Enrichment Analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; TNF, tumor necrosis factor.

GO enrichment analysis of the upregulated genes revealed significant enrichment in BP related to inflammatory response, lipopolysaccharide stimulus, response to molecules of bacterial origin, and regulation of cytokine production. The associated CC were primarily involved in cell-matrix adhesion and focal adhesions, while the MF were enriched in cytokine receptor activity and immune receptor activity (Figure 3C). KEGG pathway analysis showed that the upregulated genes were predominantly enriched in TNF signaling pathway, cytokine-cytokine receptor interaction, transcriptional misregulation in cancer, Staphylococcus aureus infection, and the HIF-1 signaling pathway (Figure 3D).

GO analysis of the downregulated genes indicated enrichment in BPs related to fatty acid metabolic processes, regulation of lipid metabolism, and monocarboxylic acid biosynthesis. MFs were enriched in lipid transporter activity and oxidoreductase activity (Figure 3E). KEGG analysis revealed that the downregulated genes were significantly enriched in adenosine triphosphate (ATP)-binding cassette (ABC) transporter signaling, adipocytokine signaling, PPAR signaling pathway, and insulin resistance (Figure 3F).

Further enrichment analysis using GSEA demonstrated that pathways related to bile acid and bile salt recycling (Figure 4A), taste transduction (Figure 4B), and ABC transporters (Figure 4C) were significantly enriched in patients with high HIF1A expression, suggesting a potential role of metabolic and transport-associated pathways in the pathogenesis of BA.

Figure 4 GSEA of the high HIF1A expression group. (A) REACTOME_RECYCLING_OF_BILE_ACIDS_AND_SALTS pathway, associated with the recycling of bile acids and salts. (B) KEGG_TASTE_TRANSDUCTION pathway, involved in taste signal transduction. (C) KEGG_ABC_TRANSPORTERS pathway, related to ABC transporters. ABC, ATP-binding cassette; ATP, adenosine triphosphate; FDR, false discovery rate; GSEA, Gene Set Enrichment Analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; NES, normalized enrichment score.

Identification of common immune-related DEGs

Patients in the validation dataset GSE15235 were classified into inflammatory and fibrotic subtypes based on molecular phenotyping. Subtype labels (inflammatory vs. fibrotic) followed the original GSE15235 annotation. This subtype information was directly obtained from the original GSE15235 dataset annotations, which defined inflammatory and fibrotic types based on histopathologic evaluation. Other clinical categories, including genetic causes (e.g., PKD1L1, ARF6, EFEMP1), CMV infection status, or cystic BA subtype, were not available in this public dataset and thus could not be analyzed. A t-test revealed that HIF1A expression was significantly higher in the inflammatory subtype compared to the fibrotic subtype (Figure 5A). Subsequently, inflammatory-type patients in GSE15235 were stratified into high and low HIF1A expression groups, resulting in the identification of 2,325 DEGs.

Figure 5 Identification and functional association analysis of common immune-related genes correlated with HIF1A. (A) In the GSE15235 dataset, HIF1A expression was significantly higher in the inflammatory subtype compared to the fibrotic subtype. (B) Venn diagram illustrating the intersection of DEGs based on HIF1A expression in the GSE46960 training set, DEGs from HIF1A-based stratification of inflammatory-type patients in GSE15235, and 782 immune infiltration-related genes, resulting in the identification of 27 common immune-related DEGs. (C) PPI network constructed via the STRING database showing extensive interactions between HIF1A (red node) and the 27 immune-related DEGs; yellow nodes represent core interacting genes. (D,E) Correlation analyses of HIF1A expression with the 27 immune-related DEGs in the GSE46960 training set (D) and the inflammatory subgroup of the GSE15235 validation set (E), indicating widespread positive associations. *, P<0.05. DEG, differentially expressed gene; PPI, protein-protein interaction.

An intersection analysis was then performed among the 461 DEGs identified from HIF1A-based grouping in the training dataset GSE46960 dataset, the 2,325 DEGs from the GSE15235 inflammatory subgroup, and a curated set of 782 immune cell infiltration-related genes. This analysis yielded 27 overlapping immune-related DEGs (Figure 5B). To ensure cross-dataset reliability, expression trends observed in the training dataset (GSE46960) were compared with those in the small independent dataset (GSE159720), which showed consistent directionality for key genes such as HIF1A. For the larger GSE15235 dataset, which includes only BA patients, analyses focused on the differences between inflammatory and fibrotic subtypes and on functional characterization after stratifying by HIF1A expression. This strategy emphasizes biological concordance and effect direction over nominal statistical significance, in line with best practices for small or control-free transcriptomic cohorts.

PPI network analysis revealed extensive potential interactions between HIF1A and the 27 immune-related DEGs (Figure 5C). Correlation analysis further confirmed that the expression of HIF1A exhibited widespread positive correlations with these 27 immune-related genes in both the training cohort and the inflammatory subgroup of the validation dataset (Figure 5D,5E).

Immune infiltration and correlation analysis

To investigate the relationship between HIF1A expression and immune cell infiltration, patients in both the training (GSE46960) and validation (GSE15235) cohorts were stratified into high and low HIF1A expression groups. Comparative analysis using the t-test revealed that, in the GSE46960 dataset, the high-HIF1A group exhibited significantly higher xCell enrichment scores for CD8⁺ T cells, CD4⁺ memory T cells, monocytes, and neutrophils compared to the low-expression group (Figure 6A). Similarly, in the inflammatory subtype of the GSE15235 dataset, higher xCell enrichment scores for macrophages and neutrophils were observed in the high HIF1A group (Figure 6B).

Figure 6 Correlation between HIF1A expression and immune cell xCell enrichment scores and immune-related gene expression. (A) In the GSE46960 dataset, the high-HIF1A group showed higher xCell enrichment scores for CD8⁺ T cells, CD4⁺ memory T cells, monocytes, and neutrophils. (B) In the inflammatory subtype of GSE15235, xCell enrichment scores were higher in the high-HIF1A group. (C,D) HIF1A expression positively correlated with neutrophil marker genes in both cohorts. (E) HIF1A expression showed significant positive correlations with the expression of multiple MHC molecules. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant (P>0.05). H, high; L, low; MHC, major histocompatibility complex; NK, natural killer.

Notably, neutrophil infiltration was consistently and significantly increased in both datasets, suggesting a potentially critical role for neutrophils in HIF1A-mediated immune responses. Further correlation analysis demonstrated that HIF1A expression was positively associated with multiple neutrophil marker genes in both the training and validation cohorts (Figure 6C,6D). In addition, HIF1A expression showed significant positive correlations with a range of major MHC molecules (Figure 6E), indicating its possible involvement in antigen presentation and adaptive immune regulation.

WGCNA

In the inflammatory subtype of the GSE15235 dataset, a total of 2,325 DEGs were identified based on HIF1A expression stratification. A weighted gene co-expression network was subsequently constructed using WGCNA, which identified several gene modules strongly associated with transcriptional traits. The network construction used a soft-threshold power (β=7) to achieve scale-free topology (R²=0.89), with a minimum module size of 10 and mergeCutHeight of 0.25, yielding stable module clustering across permutations (Figure S2). Among these, the MEblue module showed the strongest positive correlation with HIF1A expression (correlation coefficient =0.80, P=1×10⁻⁴) and was defined as the key module (Figure 7A).

Figure 7 Identification of key WGCNA co-expression modules and functional enrichment analysis based on HIF1A expression levels. (A) Heatmap showing the correlation between gene co-expression modules and HIF1A expression groups. The MEblue module exhibited the strongest positive correlation with the high HIF1A expression group (correlation coefficient =0.80, P<1×10⁻⁴), and was identified as the key module. (B) Functional enrichment analysis of genes in the MEblue module revealed significant enrichment in immune-related pathways, including cytokine signaling in the immune system (R-HSA-1280215), response to cytokine stimulus (GO:0071345), positive regulation of cytokine production (GO:0001819), and leukocyte chemotaxis (GO:0030595), among others. Different colored bars represent distinct functional categories/pathways identified by enrichment analysis. GO, Gene Ontology; WGCNA, weighted gene co-expression network analysis.

Functional enrichment analysis of genes within the MEblue module revealed significant involvement in multiple immune- and inflammation-related pathways, including cytokine signaling in the immune system, cellular responses to cytokine stimuli, positive regulation of cytokine production, leukocyte chemotaxis, and regulation of acute inflammatory responses (Figure 7B). These findings suggest that the key co-expression module associated with high HIF1A expression may contribute to the pathogenesis of BA through the regulation of immune-related BPs.

Discussion

BA is a rapidly progressive neonatal cholangiopathy with poor prognosis that often culminates in liver failure and pediatric liver transplantation. Although the precise pathogenesis remains unresolved, converging evidence implicates disordered embryogenesis, viral triggers, immune dysregulation, and genetic susceptibility. Clarifying the molecular basis of BA may enable earlier diagnosis, risk stratification, and targeted therapies (25).

Using multiple GEO cohorts and a systems framework (differential expression, immune deconvolution, pathway enrichment, and WGCNA), we prioritized candidates and modules plausibly relevant to BA biology. Across analyses, HIF1A emerged as a candidate regulator linked to inflammatory/immune programs, suggesting hypoxia-related mechanisms in BA.

DE results highlighted dysregulated genes enriched in cytokine signaling, immune activation, cell migration/adhesion, and extracellular-matrix remodeling, molecular themes concordant with clinical pathology (cholangitis, ductal obstruction, fibrosis) (9,26). The recurrent enrichment of TNF signaling and cytokine-receptor interaction pathways across datasets further validates their sustained activation and pivotal role in BA pathogenesis (8). Collectively, these molecular events likely drive the imbalance between bile duct injury and repair, progressively promoting disease progression toward end-stage liver cirrhosis (27).

Notably, this study specifically focused on the aberrant activation of development-related Hedgehog and NOTCH signaling pathways in BA, both of which play critical regulatory roles in embryonic development and bile duct morphogenesis (28,29). Integrative analysis of DEGs and developmental signaling pathways identified HIF1A as a central development-related DEG, which was significantly upregulated in BA patients. As a master transcription factor in the hypoxic response, HIF1A mediates cellular adaptive mechanisms under tissue hypoxia and regulates a variety of inflammatory cytokines, chemokines, and pro-fibrotic factors such as TGF-β, vascular endothelial growth factor (VEGF), and matrix metalloproteinases (MMPs). Its pivotal role in various hepatic fibrotic diseases has been well documented (30,31). While we used binary gene-set intersection to map DEGs to pathways, quantitative activity scoring (GSVA/ssGSEA) could further refine pathway-level inference in future analyses.

Biologically, cholestasis and ductal injury likely create a hypoxic microenvironment. In this context, HIF-1α may amplify inflammatory and fibrogenic programs, potentially accelerating progression (32).

Immune infiltration analysis further revealed that HIF1A may contribute to the pathogenesis of BA by modulating immune cell activity. Our study demonstrated a significant increase in neutrophil infiltration in BA liver tissues, with neutrophil marker gene expression positively correlating with HIF1A levels. As critical effectors of the innate immune response, neutrophils exacerbate hepatocyte injury, bile duct epithelial damage, and ECM deposition by releasing abundant inflammatory mediators, reactive oxygen species, and proteolytic enzymes, thereby playing a pivotal role in fibrosis progression (33). HIF1A may act as a key regulatory hub in neutrophil recruitment, activation, and chemotaxis, forming a hypoxia-centered inflammatory fibrotic vicious cycle (34).

Additionally, positive correlations between HIF1A and multiple major MHC molecules were observed, consistent with its involvement in adaptive immune processes and potential contribution to cellular immune activation and sustained inflammation in BA. We did not use LM22 in this study; although LM22 is blood-derived, it is widely used as an exploratory reference. Future work will triangulate the xCell findings with liver-relevant or customized signatures (e.g., CIBERSORTx, MCP-counter, EPIC) for confirmation. WGCNA identified a key gene module strongly associated with high HIF1A expression. Genes within this module were significantly enriched in leukocyte chemotaxis, cytokine signaling pathways, and immune-response regulation, indicating that BA’s molecular pathology likely involves a complex regulatory network across multiple pathways, cell types, and factors. Compared with conventional differential expression analysis, WGCNA can elucidate synergistic patterns among core regulatory factors, providing a hypothesis-generating candidate gene set for subsequent target identification and drug development (35,36). While WGCNA was performed on DEGs rather than the full transcriptome to focus on HIF1A-associated transcriptional variation, scale-free fitting and permutation testing supported network topology and module stability, lending confidence despite DEG-based prefiltering.

Moreover, this study applied FDR correction only to the differential-expression analysis, where thousands of tests were performed. For downstream analyses (enrichment, correlation, and network modeling), nominal P<0.05 was retained to balance statistical rigor and biological interpretability.

Such selective FDR application is commonly adopted in exploratory transcriptomic studies with small sample sizes, although it may increase the risk of false positives. In our study, the limited cohort size and the lack of normal controls in some GEO datasets inevitably reduce statistical power and validation depth. Moreover, the public datasets did not include structured metadata on genetic etiologies (e.g., PKD1L1, ARF6, EFEMP1), CMV infection, or cystic BA subtype, which restricted subgroup analyses by these clinically relevant categories. The small number of control samples (n=14) and the absence of uniform age or cholestasis matching further constrained between-group comparisons. We reported xCell-derived immune enrichment scores only; due to non-overlapping cell definitions between xCell and LM22/CIBERSORT, cross-method validation was not feasible without additional assumptions. Moreover, because we analyzed GEO-provided normalized matrices, array-level QC metrics (NUSE/RLE/array weights) were not recomputed; PCA was used only as an exploratory visualization. These constraints should be considered when interpreting immune-cell inferences.

Consequently, our findings should be regarded as exploratory and hypothesis-generating, pending validation in larger, clinically annotated cohorts. In addition, by the time of KPE, the extrahepatic bile ducts are already obliterated, and the sampled tissues primarily reflect secondary processes such as cholestasis, inflammation, and fibrosis rather than the initiating pathogenic mechanisms. These late-stage alterations may obscure early disease-driving molecular events, making it difficult to distinguish cause from consequence. Therefore, the present transcriptomic findings should be interpreted within this biological context, emphasizing their exploratory value and the need for future validation in early-stage clinical or experimental models. Although the sample size precluded formal multivariable or partial-correlation modeling, comparisons with independent immune indicators—including xCell-derived cell enrichment scores, neutrophil/macrophage markers consistently supported the association between HIF1A and immune activation. These results suggest that the observed relationship is not merely driven by global inflammation, but further confirmation in well-characterized cohorts is warranted.

Despite these limitations, the use of consistent preprocessing pipelines across platforms enabled exploratory comparisons and minimized obvious platform artifacts. PCA of normalized matrices was used only as an overview of sample structure.

Conclusions

In summary, this study systematically elucidated the potential bridging regulatory role of HIF1A in BA based on high-throughput transcriptomic data, and constructed a putative pathogenic model centered on HIF1A involving “hypoxia-immune activation-fibrosis”. The findings not only deepen the understanding of the immunopathological mechanisms underlying BA but also provide strong evidence supporting HIF1A and its downstream pathways as promising diagnostic biomarkers and therapeutic targets. However, certain limitations exist, including the reliance on publicly available datasets and the lack of validation using clinical samples or in vivo and in vitro functional assays. Future studies should integrate cellular and animal models to further verify the pathogenic roles of HIF1A and explore its potential as an intervention target in anti-inflammatory and anti-fibrotic therapies, thereby advancing BA management toward precision diagnosis and personalized treatment. To address the limitations of transcriptomic analyses and better distinguish primary pathogenic events from secondary inflammatory or fibrotic consequences, future studies will integrate in vitro and in vivo experimental models. Specifically, patient-derived biliary organoids and bile-duct injury animal models will be employed to test whether HIF1A activation directly promotes neutrophil recruitment, inflammatory amplification, and fibrotic remodeling. These mechanistic experiments will help validate the causal role of HIF1A and clarify whether its upregulation represents an initiating or downstream response during BA pathogenesis. In addition, clinical variables such as age at KPE, which is known to strongly influence surgical outcomes, were not available in the public datasets. Early- and late-presenting BA cases may exhibit distinct molecular and genetic features. Future studies incorporating detailed clinical metadata will be essential to determine whether age at surgery correlates with specific transcriptional profiles or disease-driving mechanisms.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the STREGA reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-590/rc

Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-590/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-590/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. All analyses were performed on de-identified, publicly available datasets from the GEO database; therefore, additional institutional review board approval and informed consent were not required.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Tam PKH, Wells RG, Tang CSM, et al. Biliary atresia. Nat Rev Dis Primers 2024;10:47. [Crossref] [PubMed]
2. Alatas FS, Lazarus G, Junaidi MC, et al. Prophylactic Antibiotics to Prevent Cholangitis in Children with Biliary Atresia After Kasai Portoenterostomy: A Meta-Analysis. J Pediatr Gastroenterol Nutr 2023;77:648-54. [Crossref] [PubMed]
3. Gopal SH, Zebda R, Mohan A, et al. Population-based screening strategies for biliary atresia in the newborn: A systematic review and meta-analysis. PLoS One 2024;19:e0307837. [Crossref] [PubMed]
4. Lendahl U, Lui VCH, Chung PHY, et al. Biliary Atresia - emerging diagnostic and therapy opportunities. EBioMedicine 2021;74:103689. [Crossref] [PubMed]
5. Fischler B, Liliemark U, Psaros Einberg A, et al. Cytomegalovirus infection in patients with biliary atresia-further questions and possible solutions. Hepatol Commun 2024;8:e0339. [Crossref] [PubMed]
6. Chen S, Li P, Wang Y, et al. Rotavirus Infection and Cytopathogenesis in Human Biliary Organoids Potentially Recapitulate Biliary Atresia Development. mBio 2020;11:e01968-20. [Crossref] [PubMed]
7. Miller PN, Baskaran S, Nijagal A. Immunology of Biliary Atresia. Semin Pediatr Surg 2024;33:151474. [Crossref] [PubMed]
8. Ye C, Zhu J, Wang J, et al. Single-cell and spatial transcriptomics reveal the fibrosis-related immune landscape of biliary atresia. Clin Transl Med 2022;12:e1070. [Crossref] [PubMed]
9. Wang J, Xu Y, Chen Z, et al. Liver Immune Profiling Reveals Pathogenesis and Therapeutics for Biliary Atresia. Cell 2020;183:1867-1883.e26. [Crossref] [PubMed]
10. Fang Z, Meng Q, Xu J, et al. Signaling pathways in cancer-associated fibroblasts: recent advances and future perspectives. Cancer Commun (Lond) 2023;43:3-41. [Crossref] [PubMed]
11. Hu Y, Bao X, Zhang Z, et al. Hepatic progenitor cell-originated ductular reaction facilitates liver fibrosis through activation of hedgehog signaling. Theranostics 2024;14:2379-95. [Crossref] [PubMed]
12. Wang T, Xu C, Zhang Z, et al. Cellular heterogeneity and transcriptomic profiles during intrahepatic cholangiocarcinoma initiation and progression. Hepatology 2022;76:1302-17. [Crossref] [PubMed]
13. Quelhas P, Breton MC, Oliveira RC, et al. HIF-1alpha-pathway activation in cholangiocytes of patients with biliary atresia: An immunohistochemical/molecular exploratory study. J Pediatr Surg 2023;58:587-94. [Crossref] [PubMed]
14. Fu M, Guo Z, Chen Y, et al. Proteomics Defines Plasma Biomarkers for the Early Diagnosis of Biliary Atresia. J Proteome Res 2024;23:1744-56. [Crossref] [PubMed]
15. Xu L, Xiao T, Zou B, et al. Identification of diagnostic biomarkers and potential therapeutic targets for biliary atresia via WGCNA and machine learning methods. Front Pediatr 2024;12:1339925. [Crossref] [PubMed]
16. Du B, Mu K, Sun M, et al. Biliary atresia and cholestasis plasma non-targeted metabolomics unravels perturbed metabolic pathways and unveils a diagnostic model for biliary atresia. Sci Rep 2024;14:15796. [Crossref] [PubMed]
17. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 2007;8:118-27. [Crossref] [PubMed]
18. Leek JT, Johnson WE, Parker HS, et al. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012;28:882-3. [Crossref] [PubMed]
19. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
20. Conesa A, Madrigal P, Tarazona S, et al. A survey of best practices for RNA-seq data analysis. Genome Biol 2016;17:13. [Crossref] [PubMed]
21. Reimand J, Isserlin R, Voisin V, et al. Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nat Protoc 2019;14:482-517. [Crossref] [PubMed]
22. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
23. Liberzon A, Birger C, Thorvaldsdóttir H, et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015;1:417-25. [Crossref] [PubMed]
24. Angel A, Naom L, Nabet-Levy S, et al. xCell 2.0: robust algorithm for cell type proportion estimation predicts response to immune checkpoint blockade. Genome Biol 2025;26:335. [Crossref] [PubMed]
25. Xiao MH, Ma D, Wu S, et al. Integrative single-cell and spatial transcriptomic analyses identify a pathogenic cholangiocyte niche and TNFRSF12A as therapeutic target for biliary atresia. Hepatology 2025;81:1146-63. [Crossref] [PubMed]
26. Li Y, Leung PS, Zhang W, et al. Immunobiology of bile and cholangiocytes. J Autoimmun 2025;151:103376. [Crossref] [PubMed]
27. de Jong IEM, Hunt ML, Chen D, et al. A fetal wound healing program after intrauterine bile duct injury may contribute to biliary atresia. J Hepatol 2023;79:1396-407. [Crossref] [PubMed]
28. Ingham PW. Hedgehog signaling. Curr Top Dev Biol 2022;149:1-58. [Crossref] [PubMed]
29. Martinez Lyons A, Boulter L. The developmental origins of Notch-driven intrahepatic bile duct disorders. Dis Model Mech 2021;14:dmm048413. [Crossref] [PubMed]
30. Chu Q, Gu X, Zheng Q, et al. Regulatory mechanism of HIF-1α and its role in liver diseases: a narrative review. Ann Transl Med 2022;10:109. [Crossref] [PubMed]
31. Luo M, Li T, Sang H. The role of hypoxia-inducible factor 1α in hepatic lipid metabolism. J Mol Med (Berl) 2023;101:487-500. [Crossref] [PubMed]
32. Zhang Y, Xu S, Shao J, et al. TRAF2 Promotes Liver Fibrosis via Regulation of the HIF-1α/GLUT1-Mediated Glycolysis in Hepatic Stellate Cells. Int J Biol Sci 2025;21:5645-65. [Crossref] [PubMed]
33. Liu K, Wang FS, Xu R. Neutrophils in liver diseases: pathogenesis and therapeutic targets. Cell Mol Immunol 2021;18:38-44. [Crossref] [PubMed]
34. Willson JA, Arienti S, Sadiku P, et al. Neutrophil HIF-1α stabilization is augmented by mitochondrial ROS produced via the glycerol 3-phosphate shuttle. Blood 2022;139:281-6. [Crossref] [PubMed]
35. Jiang W, Ye W, Tan X, et al. Network-based multi-omics integrative analysis methods in drug discovery: a systematic review. BioData Min 2025;18:27. [Crossref] [PubMed]
36. Liu Z, Lai S, Qu Q, et al. Analysis of weighted gene co-expression networks and clinical validation identify hub genes and immune cell infiltration in the endometrial cells of patients with recurrent implantation failure. Front Genet 2024;15:1292757. [Crossref] [PubMed]