# The Bioinformatic Study Uncovers Probable Critical Genes Involved in the Pathophysiology of Biliary Atresia

Authors: Si Chen; Songli Wu; Mingman Zhang
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

Background. Biliary atresia (BA) is an uncommon illness that causes the bile ducts outside and within the liver to become clogged in babies. If left untreated, the cholestasis causes increasing conjugated hyperbilirubinemia, cirrhosis, and hepatic failure. BA has a complicated aetiology, and the mechanisms that drive its development are unknown. The objective of this study was to show the role of probable critical genes involved in the pathophysiology of biliary atresia. Methods. We utilised the public Gene Expression Omnibus (GEO) microarray expression profiling dataset GSE46960 to find differentially expressed genes (DEGs) in 64 biliary atresia newborns, 14 infants with various causes of intrahepatic cholestasis, and 7 deceased-donor children as control subjects in our study. The relevant information was looked into. The important modules were identified after functional enrichment, GO and KEGG pathway analyses, protein-protein interaction (PPI) network analyses, and GSEA analysis. Results. The differential expression analysis revealed a total of 22 elevated genes. To further understand the biological activities of the DEGs, we run functional enrichment analyses on them. Meanwhile, KEGG analysis has revealed significant enrichment of pathways involved in activating cross-talking with inflammation and fibrosis in BA. SERPINE1, THBS1, CCL2, MMP7, CXCL8, EPCAM, VCAN, ITGA2, AREG, and HAS2, which may play a significant regulatory role in the pathogenesis of BA, were identified by PPI studies. Conclusion. Our findings suggested 10 hub genes and probable mechanisms of BA in the current study through bioinformatic analysis.


---

## Body

## 1. Introduction

Biliary atresia (BA) is a rare disease in which the bile ducts outside and inside the liver become blocked in newborns. Increasing evidence showed that newborn screening with direct or conjugated bilirubin results in earlier diagnosis. The serum bilirubin level after Kasai portoenterostomy is still the most accurate clinical predictor of native liver survival. Cholestasis causes increasing conjugated hyperbilirubinemia, cirrhosis, and hepatic failure if not treated [1, 2]. It is the most likely cause for a liver transplant in a youngster. BA has a tangled aetiology, with evidence pointing to viral, toxic, and genetic factors [3, 4]. The mechanisms that cause it are likewise unknown. We still do not know when BA starts or how to prevent the liver from deteriorating further [5]. A better knowledge of the aetiology of BA is required for novel therapy options other than liver transplantation to be developed. As a result, the goal of our research was to look at gene expression patterns in BA patients in order to look for potential biomarkers or pathological causes of the disease, as well as to find a better understanding and therapy for the condition.


## 2. Materials and Methods

### 2.1. Microarray Data

The gene expression profiling dataset GSE46960, which was deposited by Bessho et al. [6], was obtained using the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [7]. The dataset was created using the GPL6244 Affymetrix Human Gene 1.0 ST Array (transcript (gene) version) platform. Liver biopsy samples were obtained from 64 neonates with biliary atresia during an intraoperative cholangiogram, 14 age-matched babies with various kinds of intrahepatic cholestasis served as diseased controls, and 7 deceased-donor children served as normal controls. The age and sex of the participants, as well as their preoperative biochemical test data, were inaccessible because it was a public dataset, which looks to be a possible drawback. GPL6244’s annotation file was also obtained from the GEO.


### 2.2. Differential Expression Analysis

Using the online analytic tool GEO2R, the expression profiles of BA, non-BA patients, and healthy controls were compared to find DEGs.P values and corrected P values were calculated using T-tests. The platform’s gene probes were translated into gene names by referencing the GPL6244 platform. The genes in each sample were preserved if they matched two criteria: (1) a ∣log2fold−change∣>1 and (2) an adjusted P0.05. We identified the most important genes when the DEGs were repeated. The DEGs were found by the intersection of the two datasets, which were conducted independently for the BA versus NC and BA vs. non-BA groups.The online tool E Venn [8] (http://www.ehbio.com/test/venn/#/) was used to construct a Venn diagram of DEGs, and the heat map for the DEGs was made using the online tool xiantao Xue shu (https://www.xiantao.love/).


### 2.3. Functional Enrichment Analysis of DEGs

To improve the identification of the biological activities of DEGs, we used the web tool DAVID (https://david.ncifcrf.gov/) to conduct Gene Ontology (GO) terms [9] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [10]. Based on the GO analysis description, the gene function annotations were categorised as biological processes (BP), cellular components (CC), or molecular activities (MF). Statistical significance was defined as adjusted P values of less than 0.05. ClueGo [11, 12], a Cytoscape application (Cytoscape v3.8.0) plug-in, was utilised to demonstrate the relationship between gene enrichment analysis terms.


### 2.4. PPI Network Creation and Hub Gene Identification

A PPI network of DEGs was created using the Search Tool for the Retrieval of Interacting Genes (STRING11.5;https://string-db.org/). [13], with an interaction score cut-off of >0.4. The hub genes were found using Cytohubba [14], a Cytoscape software (Cytoscape v3.8.0) plug-in, and the key modules in the PPI network were found using molecular complex identification (MCODE 1.5.1) [15], another Cytoscape software plug-in. The DEG clustering and scoring parameters were MCODEscore=4, degreecut‐off=2, nodescorecut‐off=0.2, maxdepth=100, and k‐score=2.


### 2.5. GSEA Gene Set Enrichment Analysis

GSEA is a computer programme that determines if a set of genes that have been defined a priori demonstrate a statistically meaningful, congruent gap between different physiological situations (e.g., phenotypes). GSE46960 was submitted to Gene Set Enrichment Analysis withpermutation=1,000 using the GSEA tool (https://www.broadinstitute.org/gsea/) [16, 17]. A hypothetical P value was used to assess the statistically significant results of the enrichment score.


### 2.6. Statistical Analysis

Continuous normally distributed data are expressed as themeans±SDs. All statistical calculations were calculated through SPSS statistical software. P values < 0.05 were considered significant.


## 2.1. Microarray Data

The gene expression profiling dataset GSE46960, which was deposited by Bessho et al. [6], was obtained using the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [7]. The dataset was created using the GPL6244 Affymetrix Human Gene 1.0 ST Array (transcript (gene) version) platform. Liver biopsy samples were obtained from 64 neonates with biliary atresia during an intraoperative cholangiogram, 14 age-matched babies with various kinds of intrahepatic cholestasis served as diseased controls, and 7 deceased-donor children served as normal controls. The age and sex of the participants, as well as their preoperative biochemical test data, were inaccessible because it was a public dataset, which looks to be a possible drawback. GPL6244’s annotation file was also obtained from the GEO.


## 2.2. Differential Expression Analysis

Using the online analytic tool GEO2R, the expression profiles of BA, non-BA patients, and healthy controls were compared to find DEGs.P values and corrected P values were calculated using T-tests. The platform’s gene probes were translated into gene names by referencing the GPL6244 platform. The genes in each sample were preserved if they matched two criteria: (1) a ∣log2fold−change∣>1 and (2) an adjusted P0.05. We identified the most important genes when the DEGs were repeated. The DEGs were found by the intersection of the two datasets, which were conducted independently for the BA versus NC and BA vs. non-BA groups.The online tool E Venn [8] (http://www.ehbio.com/test/venn/#/) was used to construct a Venn diagram of DEGs, and the heat map for the DEGs was made using the online tool xiantao Xue shu (https://www.xiantao.love/).


## 2.3. Functional Enrichment Analysis of DEGs

To improve the identification of the biological activities of DEGs, we used the web tool DAVID (https://david.ncifcrf.gov/) to conduct Gene Ontology (GO) terms [9] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [10]. Based on the GO analysis description, the gene function annotations were categorised as biological processes (BP), cellular components (CC), or molecular activities (MF). Statistical significance was defined as adjusted P values of less than 0.05. ClueGo [11, 12], a Cytoscape application (Cytoscape v3.8.0) plug-in, was utilised to demonstrate the relationship between gene enrichment analysis terms.


## 2.4. PPI Network Creation and Hub Gene Identification

A PPI network of DEGs was created using the Search Tool for the Retrieval of Interacting Genes (STRING11.5;https://string-db.org/). [13], with an interaction score cut-off of >0.4. The hub genes were found using Cytohubba [14], a Cytoscape software (Cytoscape v3.8.0) plug-in, and the key modules in the PPI network were found using molecular complex identification (MCODE 1.5.1) [15], another Cytoscape software plug-in. The DEG clustering and scoring parameters were MCODEscore=4, degreecut‐off=2, nodescorecut‐off=0.2, maxdepth=100, and k‐score=2.


## 2.5. GSEA Gene Set Enrichment Analysis

GSEA is a computer programme that determines if a set of genes that have been defined a priori demonstrate a statistically meaningful, congruent gap between different physiological situations (e.g., phenotypes). GSE46960 was submitted to Gene Set Enrichment Analysis withpermutation=1,000 using the GSEA tool (https://www.broadinstitute.org/gsea/) [16, 17]. A hypothetical P value was used to assess the statistically significant results of the enrichment score.


## 2.6. Statistical Analysis

Continuous normally distributed data are expressed as themeans±SDs. All statistical calculations were calculated through SPSS statistical software. P values < 0.05 were considered significant.


## 3. Results

### 3.1. Identification of DEGs in Biliary Atresia

The gene expression profile of the GSE46960 dataset comprised data from three separate groups (Table1). Using a fold change (FC) value of ∣log2FC∣>2 and a P value of 0.05 as the cut-off, a total of 22 DEGs, all upregulated genes, were obtained from the notably regulated gene in the biliary atresia samples over both diseased control (non-BA) and normal control (NC) samples (Table 2). For the distribution of DEGs, an online tool was utilised to construct a Venn diagram and heat maps (Figures 1(a) and 1(b)).Table 1
Details GEO biliary atresia data.
SampleGEOPlatformBANon-BANCLiverGSE46960GPL624464147Table 2
List of genes specifically regulated in the biliary atresia samples over both diseased control (non-BA) and normal control (NC).
Gene symbolAdjustedP valueFold changeDescriptionBA-NCBA-non-BABA-NCBA-non-BAUpregulated genesEPCAM4.73E-254.56E-054.371.29Epithelial cell adhesion moleculeSPP18.86E-174.56E-053.261.37Secreted phosphoprotein 1ANKRD14.14E-101.73E-063.112.05Ankyrin repeat domain 1MMP73.76E-088.21E-052.941.85Matrix metallopeptidase 7LUM3.21E-126.17E-052.701.26LumicanRGS47.25E-125.86E-042.541.03Regulator of G-protein signaling 4EMP11.16E-051.96E-032.401.41Epithelial membrane protein 1CFTR1.41E-082.82E-032.331.05Cystic fibrosis transmembrane conductance regulatorHAS23.69E-061.62E-042.281.56Hyaluronan synthase 2KRT231.34E-085.29E-042.191.14Keratin 23VCAN6.40E-101.52E-052.121.32VersicanCXCL83.55E-094.56E-051.931.20C-X-C motif chemokine ligand 8CCL201.75E-031.14E-021.871.24C-C motif chemokine ligand 20VTCN14.14E-081.11E-041.811.09V-set domain containing T cell activation inhibitor 1ITGA26.29E-071.50E-041.681.10Integrin subunit alpha 2SERPINE11.28E-029.47E-031.531.36Serpin family E member 1THBS15.61E-063.54E-041.521.00Thrombospondin 1CCL22.51E-044.05E-051.381.37C-C motif chemokine ligand 2TM4SF11.14E-031.52E-051.381.71Transmembrane 4 L six family member 1LAMC21.99E-064.56E-051.351.02Laminin subunit gamma 2AREG6.89E-032.44E-031.341.28AmphiregulinSLC2A31.05E-022.66E-031.151.16Solute carrier family 2 member 3Figure 1
(a) The selection of 22 genes based on BA vs. NC and non-BA as sick controls. (b) The expression data is displayed as a data matrix, with each row corresponding to a gene and each column to a sample. The colour ratio of the upper left corner is used to convey the amount of emotion. The top tree view demonstrates hierarchical clustering and indicates the degree of gene expression relatedness. Abbreviations: DEG: differentially expressed genes; BA: biliary atresia; non-BA: other causes of intrahepatic cholestasis except biliary atresia; NC: normal control; FC: fold change.
(a)(b)


### 3.2. GO and KEGG Pathway Analysis for Identifying the DEGs

DEGs were studied using the DAVID online tool for functional and pathway enrichment. The most three important processes revealed by GO analysis among the annotations of BP were extracellular matrix organization, cellular response to tumor necrosis factor, and cell adhesion. The three most important processes revealed among the CC annotations were extracellular space, extracellular area, and extracellular exosome. Finally, the three most significant processes among the MF annotations were receptor binding, glycosaminoglycan binding, and cytokine activity. Table3 and Figure 2 show the number of genes and P values of the top 8 enriched functional words based on the criteria.Table 3
Top 8 of the most significantly enriched GO terms.
Pathway IDTermsGene countP valueBPGO:0030198Extracellular matrix organization75.60E-08GO:0071356Cellular response to tumor necrosis factor61.19E-07GO:0071347Cellular response to interleukin-151.08E-06GO:0071222Cellular response to lipopolysaccharide56.97E-06GO:0007155Cell adhesion61.26E-04GO:0006954Inflammatory response57.57E-04GO:0008284Positive regulation of cell proliferation51.63E-03GO:0006955Immune response41.12E-02CCGO:0005615Extracellular space117.80E-07GO:0005576Extracellular region103.53E-05GO:0009986Cell surface63.13E-04GO:0031012Extracellular matrix53.28E-04GO:0009897External side of plasma membrane32.47E-02GO:0070062Extracellular exosome83.29E-02GO:0048471Perinuclear region of cytoplasm43.32E-02GO:0005578Proteinaceous extracellular matrix33.77E-02MFGO:0005539Glycosaminoglycan binding59.29E-05GO:0005125Cytokine activity51.20E-04GO:0005102Receptor binding91.51E-04GO:0008201Heparin binding48.46E-04GO:0008009Chemokine activity31.52E-03GO:1901681Sulfur compound binding42.49E-03GO:0043236Laminin-binding22.92E-02GO:0050840Extracellular matrix binding23.04E-02Notes. GO: Gene Ontology; BP: biological progress; CC: cellular component; MF: molecular function.Figure 2
Results of GO enrichment. The ordinate shows the number and ratio of differentially expressed genes, whereas the abscissa reflects the enriched GO. Biological process, cellular component, and molecular function are all represented by distinct colours. Abbreviation: GO: gene ontology.The DEGs’ cell signaling pathway enrichment study yielded a total of eight relevant pathways that were investigated. ECM-receptor interaction, malaria, the PI3K-Akt signaling pathway, and others were among the cellular signaling pathways linked to biliary atresia. Table4 and Figure 3 describe the specific enriched pathways discovered by DEG analysis. Figure 4 shows the relationship between the words of gene enrichment analysis.Table 4
Significantly enriched KEGG pathway.
Pathway IDTermGene countP valueGeneshsa04512ECM-receptor interaction45.10E-04THBS1, ITGA2, LAMC2, SPP1hsa05144Malaria33.69E-03CCL2, THBS1, CXCL8hsa04510Focal adhesion46.06E-03THBS1, ITGA2, LAMC2, SPP1hsa05323Rheumatoid arthritis31.15E-02CCL20, CCL2, CXCL8hsa05142Chagas disease (American trypanosomiasis)31.58E-02CCL2, CXCL8, SERPINE1hsa04151PI3K-Akt signaling pathway42.46E-02THBS1, ITGA2, LAMC2, SPP1hsa04062Chemokine signaling pathway34.66E-02CCL20, CCL2, CXCL8Notes. KEGG: Kyoto Encyclopedia of Genes and Genomes.Figure 3
KEGG pathway analysis of the differentially expressed genes in BA. Abbreviations: KEGG: Kyoto Encyclopedia of Genes and Genomes; BA: biliary atresia.Figure 4
Correlation between terms of gene enrichment analysis.


### 3.3. Construction of the PPI Network and Identification of Hub Genes

STRING (https://string-db.org/) is a public database that contains known and predicted protein interactions. PPI [18] is important for studying protein function since it can help elucidate the role of protein control. STRING’s official website was used to submit the 22 DEGs from the GSE46960 dataset in order to get protein interrelationships. The minimum required interaction score was set at 0.15 in order to see the interaction networks with Cytoscape (version v3.9.0) [19]. There were 22 nodes and 107 edges in the PPI network that resulted. The network visualisation created using STRING’s official website is shown in Figure 5(a). The degree of linkage between DEGs and genes was used to screen for hub genes, and the DEGs with the ten highest degrees were identified as hub genes (Table 5 and Figure 5(b)).Figure 5
The PPI network and the most significant modules of DEGs. (a) String software was used to examine the PPI network. In the PPI network, there were 22 nodes and 120 edges. (b) CytoHubba found the most important module. Abbreviation: DEG: differentially expressed gene; PPI: protein-protein interaction.
(a)(b)Table 5
Top 10 hub genes ranked by Cytobubba.
Gene symbolScoreFull nameDistributionSERPINE132Serpin family E member 1Serine protease inhibitor; plasminogen activator inhibitor 1. Tissue plasminogen activator, urokinase, protein C, and matriptase-3/TMPRSS7 are all attracted to this inhibitor. Its quick interaction with PLAT might be a crucial control point in fibrinolysis regulation.THBS132Thrombospondin 1Cell-to-cell and cell-to-matrix adhesion glycoprotein thrombospondin-1 modulates cell-to-cell and cell-to-matrix connections. This chemical binds to heparin. Dentin and dental pulp dentinogenesis and/or maintenance may be affected (by similarity). A CD36 ligand is responsible for the antiangiogenic actions. It has a function in the ER stress response by interacting with the activating transcription factor 6 alpha (ATF6), which produces adaptive ER stress response factors (by similarity)CCL230C-C motif chemokine ligand 2Chemotactic factor with a C-C motif, which attracts monocytes and basophils but not neutrophils or eosinophils. Increases the antitumor activity of monocytes. Has been linked to the development of disorders with monocytic infiltrates, such as psoriasis, rheumatoid arthritis, and atherosclerosis. During the disease process of atherosclerosis, it may be implicated in the recruitment of monocytes into the artery wall.MMP730Matrix metallopeptidase 7Matrilysin degrades casein, gelatins of types I, III, IV, and V, and fibronectin. Activates procollagenase; M10 matrix metallopeptidasesCXCL830C-X-C motif chemokine ligand 8IL-8 is a chemotactic factor that attracts neutrophils, basophils, and T cells but not monocytes. It has a role in neutrophil activation as well. In response to an inflammatory stimulation, it is secreted by a variety of cell types. When compared to IL-8 (1-77), IL-8 (6-77) has a 5-10-fold higher activity on neutrophil activation, IL-8 (5-77) has enhanced activity on neutrophil activation, and IL-8 (7-77) has a stronger affinity to CXCR1 and CXCR2 receptors.EPCAM28Epithelial cell adhesion moleculeIntestinal epithelial cells (IECs) and intraepithelial lymphocytes (IELs) establish a physical homophilic contact molecule at the mucosal epithelium to build an immunological barrier as a first line of defence against mucosal infection. Plays a function in the proliferation and differentiation of embryonic stem cells. The expression of FABP5, MYC, and cyclins A and E is all increased.VCAN24VersicanIntercellular signaling and cell-extracellular matrix interactions may be mediated by the Versican core protein. It might play a part in cell motility, proliferation, and differentiation. Has a C-type lectin domain that binds hyaluronic acid.ITGA224Integrin subunit alpha 2Integrin alpha-2 (integrin alpha-2/beta-1) is a receptor for laminin, collagen, collagen C-propeptides, fibronectin, and E-cadherin. It recognises the proline-hydroxylated sequence G-F-P-G-E-R in collagen. Platelet and other cell adhesion to collagens, collagen and collagenase gene expression control, force creation, and structuring of newly formed extracellular matrix are all mediated by CD molecules.AREG24AmphiregulinAmphiregulin is a ligand for the EGF receptor/EGFR. Amphiregulin is an amitogen and autocrine growth factor for a range of target cells, including astrocytes, Schwann cells, and fibroblasts.HAS222Hyaluronan synthase 2Hyaluronan synthase 2 adds GlcNAc or GlcUA monosaccharides to the nascent hyaluronan polymer. As a result, it is essential for the formation of hyaluronan, a major component of most extracellular matrices that regulates cell adhesion, migration, and differentiation and plays a structural role in tissue architecture. This is one of the isozymes that catalyses the process and is responsible for the production of high-molecular-mass hyaluronan. A key phase in the creation of the heart is the conversion of endocardial cushion cells to mesenchymal cells.


### 3.4. GSEA Analysis of All Detected Genes

GSEA was used to find gene sets with a statistically significant difference between BA and NC participants, and it revealed that the BA subjects had the highest enriched gene sets of all discovered genes. There are 4662 gene symbols in the dataset, 3651 of which are elevated in phenotype BA. ECM receptor interaction, integrin cell-surface interactions, and Andersen cholangiocarcinoma class1 were the top six most significant-enriched gene sets positively correlated with the BA subjects, followed by uterine fibroid up, ECM proteoglycans, and nonintegrin membrane ECM interactions (Figures6(a)–6(f)).Figure 6
GSEA plot showing the most enriched gene sets of all detected genes in the BA subjects. ECM receptor interaction (a), integrin cell surface contacts (b), and cholangiocarcinoma class1 are the top-six most significantly upregulated enriched gene sets in BA individuals (c). ECM proteoglycans (d), uterine fibroid uptake (e), nonintegrin membrane ECM interactions (f). Abbreviation: GSEA: gene set enrichment analysis; NES: normalized enrichment score.
(a)(b)(c)(d)(e)(f)


## 3.1. Identification of DEGs in Biliary Atresia

The gene expression profile of the GSE46960 dataset comprised data from three separate groups (Table1). Using a fold change (FC) value of ∣log2FC∣>2 and a P value of 0.05 as the cut-off, a total of 22 DEGs, all upregulated genes, were obtained from the notably regulated gene in the biliary atresia samples over both diseased control (non-BA) and normal control (NC) samples (Table 2). For the distribution of DEGs, an online tool was utilised to construct a Venn diagram and heat maps (Figures 1(a) and 1(b)).Table 1
Details GEO biliary atresia data.
SampleGEOPlatformBANon-BANCLiverGSE46960GPL624464147Table 2
List of genes specifically regulated in the biliary atresia samples over both diseased control (non-BA) and normal control (NC).
Gene symbolAdjustedP valueFold changeDescriptionBA-NCBA-non-BABA-NCBA-non-BAUpregulated genesEPCAM4.73E-254.56E-054.371.29Epithelial cell adhesion moleculeSPP18.86E-174.56E-053.261.37Secreted phosphoprotein 1ANKRD14.14E-101.73E-063.112.05Ankyrin repeat domain 1MMP73.76E-088.21E-052.941.85Matrix metallopeptidase 7LUM3.21E-126.17E-052.701.26LumicanRGS47.25E-125.86E-042.541.03Regulator of G-protein signaling 4EMP11.16E-051.96E-032.401.41Epithelial membrane protein 1CFTR1.41E-082.82E-032.331.05Cystic fibrosis transmembrane conductance regulatorHAS23.69E-061.62E-042.281.56Hyaluronan synthase 2KRT231.34E-085.29E-042.191.14Keratin 23VCAN6.40E-101.52E-052.121.32VersicanCXCL83.55E-094.56E-051.931.20C-X-C motif chemokine ligand 8CCL201.75E-031.14E-021.871.24C-C motif chemokine ligand 20VTCN14.14E-081.11E-041.811.09V-set domain containing T cell activation inhibitor 1ITGA26.29E-071.50E-041.681.10Integrin subunit alpha 2SERPINE11.28E-029.47E-031.531.36Serpin family E member 1THBS15.61E-063.54E-041.521.00Thrombospondin 1CCL22.51E-044.05E-051.381.37C-C motif chemokine ligand 2TM4SF11.14E-031.52E-051.381.71Transmembrane 4 L six family member 1LAMC21.99E-064.56E-051.351.02Laminin subunit gamma 2AREG6.89E-032.44E-031.341.28AmphiregulinSLC2A31.05E-022.66E-031.151.16Solute carrier family 2 member 3Figure 1
(a) The selection of 22 genes based on BA vs. NC and non-BA as sick controls. (b) The expression data is displayed as a data matrix, with each row corresponding to a gene and each column to a sample. The colour ratio of the upper left corner is used to convey the amount of emotion. The top tree view demonstrates hierarchical clustering and indicates the degree of gene expression relatedness. Abbreviations: DEG: differentially expressed genes; BA: biliary atresia; non-BA: other causes of intrahepatic cholestasis except biliary atresia; NC: normal control; FC: fold change.
(a)(b)


## 3.2. GO and KEGG Pathway Analysis for Identifying the DEGs

DEGs were studied using the DAVID online tool for functional and pathway enrichment. The most three important processes revealed by GO analysis among the annotations of BP were extracellular matrix organization, cellular response to tumor necrosis factor, and cell adhesion. The three most important processes revealed among the CC annotations were extracellular space, extracellular area, and extracellular exosome. Finally, the three most significant processes among the MF annotations were receptor binding, glycosaminoglycan binding, and cytokine activity. Table3 and Figure 2 show the number of genes and P values of the top 8 enriched functional words based on the criteria.Table 3
Top 8 of the most significantly enriched GO terms.
Pathway IDTermsGene countP valueBPGO:0030198Extracellular matrix organization75.60E-08GO:0071356Cellular response to tumor necrosis factor61.19E-07GO:0071347Cellular response to interleukin-151.08E-06GO:0071222Cellular response to lipopolysaccharide56.97E-06GO:0007155Cell adhesion61.26E-04GO:0006954Inflammatory response57.57E-04GO:0008284Positive regulation of cell proliferation51.63E-03GO:0006955Immune response41.12E-02CCGO:0005615Extracellular space117.80E-07GO:0005576Extracellular region103.53E-05GO:0009986Cell surface63.13E-04GO:0031012Extracellular matrix53.28E-04GO:0009897External side of plasma membrane32.47E-02GO:0070062Extracellular exosome83.29E-02GO:0048471Perinuclear region of cytoplasm43.32E-02GO:0005578Proteinaceous extracellular matrix33.77E-02MFGO:0005539Glycosaminoglycan binding59.29E-05GO:0005125Cytokine activity51.20E-04GO:0005102Receptor binding91.51E-04GO:0008201Heparin binding48.46E-04GO:0008009Chemokine activity31.52E-03GO:1901681Sulfur compound binding42.49E-03GO:0043236Laminin-binding22.92E-02GO:0050840Extracellular matrix binding23.04E-02Notes. GO: Gene Ontology; BP: biological progress; CC: cellular component; MF: molecular function.Figure 2
Results of GO enrichment. The ordinate shows the number and ratio of differentially expressed genes, whereas the abscissa reflects the enriched GO. Biological process, cellular component, and molecular function are all represented by distinct colours. Abbreviation: GO: gene ontology.The DEGs’ cell signaling pathway enrichment study yielded a total of eight relevant pathways that were investigated. ECM-receptor interaction, malaria, the PI3K-Akt signaling pathway, and others were among the cellular signaling pathways linked to biliary atresia. Table4 and Figure 3 describe the specific enriched pathways discovered by DEG analysis. Figure 4 shows the relationship between the words of gene enrichment analysis.Table 4
Significantly enriched KEGG pathway.
Pathway IDTermGene countP valueGeneshsa04512ECM-receptor interaction45.10E-04THBS1, ITGA2, LAMC2, SPP1hsa05144Malaria33.69E-03CCL2, THBS1, CXCL8hsa04510Focal adhesion46.06E-03THBS1, ITGA2, LAMC2, SPP1hsa05323Rheumatoid arthritis31.15E-02CCL20, CCL2, CXCL8hsa05142Chagas disease (American trypanosomiasis)31.58E-02CCL2, CXCL8, SERPINE1hsa04151PI3K-Akt signaling pathway42.46E-02THBS1, ITGA2, LAMC2, SPP1hsa04062Chemokine signaling pathway34.66E-02CCL20, CCL2, CXCL8Notes. KEGG: Kyoto Encyclopedia of Genes and Genomes.Figure 3
KEGG pathway analysis of the differentially expressed genes in BA. Abbreviations: KEGG: Kyoto Encyclopedia of Genes and Genomes; BA: biliary atresia.Figure 4
Correlation between terms of gene enrichment analysis.


## 3.3. Construction of the PPI Network and Identification of Hub Genes

STRING (https://string-db.org/) is a public database that contains known and predicted protein interactions. PPI [18] is important for studying protein function since it can help elucidate the role of protein control. STRING’s official website was used to submit the 22 DEGs from the GSE46960 dataset in order to get protein interrelationships. The minimum required interaction score was set at 0.15 in order to see the interaction networks with Cytoscape (version v3.9.0) [19]. There were 22 nodes and 107 edges in the PPI network that resulted. The network visualisation created using STRING’s official website is shown in Figure 5(a). The degree of linkage between DEGs and genes was used to screen for hub genes, and the DEGs with the ten highest degrees were identified as hub genes (Table 5 and Figure 5(b)).Figure 5
The PPI network and the most significant modules of DEGs. (a) String software was used to examine the PPI network. In the PPI network, there were 22 nodes and 120 edges. (b) CytoHubba found the most important module. Abbreviation: DEG: differentially expressed gene; PPI: protein-protein interaction.
(a)(b)Table 5
Top 10 hub genes ranked by Cytobubba.
Gene symbolScoreFull nameDistributionSERPINE132Serpin family E member 1Serine protease inhibitor; plasminogen activator inhibitor 1. Tissue plasminogen activator, urokinase, protein C, and matriptase-3/TMPRSS7 are all attracted to this inhibitor. Its quick interaction with PLAT might be a crucial control point in fibrinolysis regulation.THBS132Thrombospondin 1Cell-to-cell and cell-to-matrix adhesion glycoprotein thrombospondin-1 modulates cell-to-cell and cell-to-matrix connections. This chemical binds to heparin. Dentin and dental pulp dentinogenesis and/or maintenance may be affected (by similarity). A CD36 ligand is responsible for the antiangiogenic actions. It has a function in the ER stress response by interacting with the activating transcription factor 6 alpha (ATF6), which produces adaptive ER stress response factors (by similarity)CCL230C-C motif chemokine ligand 2Chemotactic factor with a C-C motif, which attracts monocytes and basophils but not neutrophils or eosinophils. Increases the antitumor activity of monocytes. Has been linked to the development of disorders with monocytic infiltrates, such as psoriasis, rheumatoid arthritis, and atherosclerosis. During the disease process of atherosclerosis, it may be implicated in the recruitment of monocytes into the artery wall.MMP730Matrix metallopeptidase 7Matrilysin degrades casein, gelatins of types I, III, IV, and V, and fibronectin. Activates procollagenase; M10 matrix metallopeptidasesCXCL830C-X-C motif chemokine ligand 8IL-8 is a chemotactic factor that attracts neutrophils, basophils, and T cells but not monocytes. It has a role in neutrophil activation as well. In response to an inflammatory stimulation, it is secreted by a variety of cell types. When compared to IL-8 (1-77), IL-8 (6-77) has a 5-10-fold higher activity on neutrophil activation, IL-8 (5-77) has enhanced activity on neutrophil activation, and IL-8 (7-77) has a stronger affinity to CXCR1 and CXCR2 receptors.EPCAM28Epithelial cell adhesion moleculeIntestinal epithelial cells (IECs) and intraepithelial lymphocytes (IELs) establish a physical homophilic contact molecule at the mucosal epithelium to build an immunological barrier as a first line of defence against mucosal infection. Plays a function in the proliferation and differentiation of embryonic stem cells. The expression of FABP5, MYC, and cyclins A and E is all increased.VCAN24VersicanIntercellular signaling and cell-extracellular matrix interactions may be mediated by the Versican core protein. It might play a part in cell motility, proliferation, and differentiation. Has a C-type lectin domain that binds hyaluronic acid.ITGA224Integrin subunit alpha 2Integrin alpha-2 (integrin alpha-2/beta-1) is a receptor for laminin, collagen, collagen C-propeptides, fibronectin, and E-cadherin. It recognises the proline-hydroxylated sequence G-F-P-G-E-R in collagen. Platelet and other cell adhesion to collagens, collagen and collagenase gene expression control, force creation, and structuring of newly formed extracellular matrix are all mediated by CD molecules.AREG24AmphiregulinAmphiregulin is a ligand for the EGF receptor/EGFR. Amphiregulin is an amitogen and autocrine growth factor for a range of target cells, including astrocytes, Schwann cells, and fibroblasts.HAS222Hyaluronan synthase 2Hyaluronan synthase 2 adds GlcNAc or GlcUA monosaccharides to the nascent hyaluronan polymer. As a result, it is essential for the formation of hyaluronan, a major component of most extracellular matrices that regulates cell adhesion, migration, and differentiation and plays a structural role in tissue architecture. This is one of the isozymes that catalyses the process and is responsible for the production of high-molecular-mass hyaluronan. A key phase in the creation of the heart is the conversion of endocardial cushion cells to mesenchymal cells.


## 3.4. GSEA Analysis of All Detected Genes

GSEA was used to find gene sets with a statistically significant difference between BA and NC participants, and it revealed that the BA subjects had the highest enriched gene sets of all discovered genes. There are 4662 gene symbols in the dataset, 3651 of which are elevated in phenotype BA. ECM receptor interaction, integrin cell-surface interactions, and Andersen cholangiocarcinoma class1 were the top six most significant-enriched gene sets positively correlated with the BA subjects, followed by uterine fibroid up, ECM proteoglycans, and nonintegrin membrane ECM interactions (Figures6(a)–6(f)).Figure 6
GSEA plot showing the most enriched gene sets of all detected genes in the BA subjects. ECM receptor interaction (a), integrin cell surface contacts (b), and cholangiocarcinoma class1 are the top-six most significantly upregulated enriched gene sets in BA individuals (c). ECM proteoglycans (d), uterine fibroid uptake (e), nonintegrin membrane ECM interactions (f). Abbreviation: GSEA: gene set enrichment analysis; NES: normalized enrichment score.
(a)(b)(c)(d)(e)(f)


## 4. Discussion

Biliary atresia (BA) is a fibroinflammatory disease of the intra- and extrahepatic biliary tree. In order to have a better understanding of the underlying cause(s) and pathogenesis of the disease, the National Institutes of Diabetes and Digestive and Kidney Diseases sponsored researches that study the promising and innovative approaches. In this investigation, we used the GEO database to screen for DEGs and acquire gene expression profiles from patients with BA, non-BA, and normal controls. There were a total of 22 DEGs confirmed.The DEGs were considerably enriched in the cellular response to interleukin-1, according to BP in GO annotation, which was consistent with earlier evidence that inhibiting IL-1-mediated inflammation may be advantageous in selective liver fibrotic disease [20]. Other enhanced gene sets of DEGs in the BP of GO, such as immunological and inflammatory responses, have been linked to biliary atresia [21, 22]. The extracellular exosome was shown to be rich in CC.Exosomes have been explored as disease biomarkers [23, 24] or cell-cell communication factors because of their role in carrying a variety of proteins, noncoding RNA, and coding RNA from different cells. A new study suggests that serum exosomal H19 might be exploited as a noninvasive diagnostic biomarker and treatment target for BA [25]. According to KEGG enrichment analysis, DEGs are also detected in the ECM-receptor interaction, focal adhesion, PI3K-Akt signaling pathway, and chemokine signaling pathway. All of these results corroborated previous findings that BA interacts with inflammation and fibrosis [26].SERPINE1, THBS1, CCL2, MMP7, CXCL8, EPCAM, VCAN, ITGA2, AREG, and HAS2 were among the 10 hub genes discovered in this study. Interleukin- (IL-) 8 (CXCL8) may mediate liver damage in BA by enhancing ductular response and related hepatic fibrogenesis, according to Godbole et al. [27]. The serum MMP-7 test, according to Yang et al., shows excellent sensitivity and specificity for distinguishing BA from other newborn cholestasis and may be a valid biomarker for BA [28]. SERPINE1 can be targeted to prevent biliary fibrosis, according to Aseem et al.The most significant-enriched gene set connected with the BA individuals, according to GSEA, was ECM receptor interaction. Many studies have linked oxidative stress to liver fibrosis. It has been discovered that ROS can activate KCs (Kupffer cells) to trigger the inflammatory response, which subsequently leads to HSC (activated hepatic stellate cells) activation to create ECM proteins [29, 30] and fibrosis. It will offer a fresh look at the treatment strategy for BA’s fibrosis mechanism. The limit of this study is that there are only bioinformatic analysis and did not have cell and animal experiments. Therefore, many investigations need to be added to the article.


## 5. Conclusion

With bioinformatic analysis, we found 10 hub genes and probable mechanisms of BA in the current study. More research is needed to confirm the hub genes and identify relevant processes. All of the findings will pave the way for a possible treatment strategy for biliary atresia and associated fibrotic illnesses.



---

*Source: 9108804-2022-06-21.xml*