Functional investigation and two-sample Mendelian randomization study of primary biliary cholangitis hub genes

Abstract

BACKGROUND

The identification of specific gene expression patterns is crucial for understanding the mechanisms underlying primary biliary cholangitis (PBC) and finding relevant biomarkers for diagnosis and therapeutic evaluation.

AIM

To determine PBC-associated hub genes and assess their clinical utility for disease prediction.

METHODS

PBC expression data were obtained from the Gene Expression Omnibus database. Overlapping genes from differential expression analysis and weighted gene co-expression network analysis (WGCNA) were identified as key genes for PBC. Kyoto Encyclopedia of Genes and Genomes and Gene Ontology analyses were performed to explore the potential roles of key genes. Hub genes were identified in protein-protein interaction (PPI) networks using the Degree algorithm in Cytoscape software. The relationship between hub genes and immune cells was investigated. Finally, a Mendelian randomization study was conducted to determine the causal effects of hub genes on PBC.

RESULTS

We identified 71 overlapping key genes using differential expression analysis and WGCNA. These genes were primarily enriched in pathways related to cytokine-cytokine receptor interaction, and Th1, Th2, and Th17 cell differentiation. We utilized Cytoscape software and identified five hub genes (CD247, IL10, CCL5, CCL3, and STAT3) in PPI networks. These hub genes showed a strong correlation with immune cell infiltration in PBC. However, inverse variance weighting analysis did not indicate the causal effects of hub genes on PBC risk.

CONCLUSION

Hub genes can potentially serve as valuable biomarkers for PBC prediction and treatment, thereby offering significant clinical utility.

Key Words: Primary biliary cholangitis; Weighted gene co-expression network analysis; Hub genes; Mendelian randomization; Bioinformatic analysis

Core Tip: This study identified five hub genes (CD247, IL10, CCL5, CCL3, and STAT3) associated with primary biliary cholangitis (PBC) through comprehensive bioinformatics analysis. These hub genes were enriched in immune-related pathways and strongly correlated with immune cell infiltration in PBC. Although hub genes did not have a causal effect on PBC risk, they provided valuable insights into the molecular mechanisms of PBC and showed potential as biomarkers for PBC prediction and treatment.

INTRODUCTION

Primary biliary cholangitis (PBC) is a chronic autoimmune liver disease distinguished by the gradual destruction of small bile ducts within the liver, leading to impaired bile flow, toxic bile acid accumulation, and liver damage[1,2]. Several challenges are associated with the medical management of PBC. Diagnosing PBC can be challenging, as early-stage disease may be asymptomatic or occur with nonspecific symptoms, resulting in delayed diagnosis and treatment initiation. Many patients do not respond to ursodeoxycholic acid therapy, and alternative treatment options are limited. Finally, the long-term prognosis of PBC remains uncertain, as the process of disease progression can vary widely among individuals[3]. Despite the advancements in the understanding the disease and availability of treatment options, several challenges are associated with effective PBC diagnosis and management. Further research is needed to improve early detection, develop alternative therapies, and enhance overall disease management[4].

A crucial aspect in unraveling the intricate microscopic mechanisms underlying a wide range of disorders, particularly PBC, lies in the precise characterization of specific gene expression patterns. This process is vital for identifying crucial biomarkers that facilitate diagnostic accuracy and therapeutic evaluation[5]. Sophisticated bioinformatics techniques, including weighted gene co-expression network analysis (WGCNA), Gene Ontology (GO) annotations, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, and gene set enrichment analysis have been tailored to decipher complex disease pathways. WGCNA serves as a valuable tool for precisely identifying therapeutic targets and potential biomarkers[6,7]. Hence, the core purpose of this investigation is to elucidate novel genetic entities, biomarkers, and the fundamental mechanisms underlying PBC.

In recent years, Mendelian randomization (MR) has been firmly established as a potent tool for investigating potential causal relationships. MR uses single nucleotide polymorphisms (SNPs) as instrumental variables (IVs) to leverage the natural randomization of genetic variants and evaluate the causal effects of exposure factors on outcomes[8]. This methodology boasts numerous advantages, including the ability to produce robust causal inferences, effectively controlling for confounding factors, and mitigating the bias resulting from reverse causality[9].

In this study, we analyzed liver tissue samples from patients diagnosed with PBC and compared them to samples from healthy controls (HCs). Our primary objective was to identify genes exhibiting differential expression patterns and examine their relationship with PBC. Using WGCNA, we identified the most relevant PBC-associated gene modules, enabling us to select genes for further investigation. Next, we identified five key genes-CD247, IL10, CCL5, CCL3, and STAT3-as promising diagnostic biomarkers for PBC from the protein-protein interaction (PPI) network. These genes could significantly improve diagnostic capabilities and elucidate the mechanisms underlying PBC. Additionally, we conducted MR analyses to investigate potential causal relationships between these five hub genes and PBC.

MATERIALS AND METHODS

Data source

The GSE79850 dataset was accessed and downloaded from the Gene Expression Omnibus (GEO) database. These datasets consisted of clinical features and gene expression data obtained from liver tissue samples of individuals diagnosed with PBC and HC. GSE79850 was comprised of 8 HC samples and 16 PBC patient samples. This dataset was generated through transcriptome microarray analysis of liver tissues obtained from individuals with PBC and HCs.

Identification of differentially expressed genes

Initially, we utilized R software (version 4.3.0) to read data from the GSE79850 dataset and performed preprocessing steps, including normalization. Subsequently, we conducted differentially expressed gene (DEG) analyses using the “limma” package, to identify genes with significant differences between the PBC and HC groups. For testing multiple corrections, we adopted the Benjamini-Hochberg method, a widely accepted approach for controlling the false discovery rate. The Benjamini-Hochberg method adjusts the raw P values to account for the increased likelihood of false positives while simultaneously performing many statistical tests. Specifically, we defined a gene as being differentially expressed if it exhibited a |logFC| ≥ 1 and an adjusted P value ≤ 0.05. After assessing the significance of expression levels, volcano plots, and DEG expression heatmaps were generated using the R packages “ggplot2” and "pheatmap", respectively.

Weighted gene co-expression network analysis

WGCNA is a systematic biological approach commonly used to characterize genetic association patterns across samples, aiming to identify highly interconnected gene modules and their relationship with phenotypes. It has been successfully used to identify candidate markers. In our study, we utilized the “WGCNA” R package to construct a gene co-expression network specific to PBC[6]. Subsequently, we assessed the correlation between different modules and the pathogenic mechanism of PBC, ultimately selecting the most relevant module as central genes derived from WGCNA.

GO/KEGG analysis

The overlapping genes from DEG analysis and WGCNA were identified as key genes in PBC pathogenesis. Subsequently, we conducted GO and KEGG enrichment analysis using the “clusterProfiler” R package to elucidate the mechanisms underlying disease progression and pathogenesis in PBC[10].

Discovery of hub genes in the PPI network

To explore the PPI networks and gene interactions relevant to PBC, we utilized the STRING database. The Degree algorithm in Cytoscape software facilitated the evaluation and ranking of the importance of genes within these networks. In a PPI network, the degree value reflects the number of interactions a gene has with other genes in the network, thus serving as an indicator of its significance within the network. This approach allowed us to identify hub genes critical to the gene interactions associated with the disease.

Nomogram model construction

To construct a nomogram model for predicting the PBC risk, we utilized the “rms” package[11]. The predictive power of the nomogram model was evaluated using Harrell’s concordance index. We employed the “ROC” package to construct a receiver operator characteristic (ROC) curve and assess the diagnostic efficacy of candidate biomarkers. The accuracy of the model was determined by calculating the area under the ROC curve (AUC), with 0.9 ≤ AUC < 1 indicating excellent accuracy.

Immune cell analysis

In order to investigate the role of immune cells in PBC, we conducted a CIBERSORT analysis to compare the immune cell infiltration levels of 22 immune cell types between the PBC and HC groups[12]. Furthermore, we evaluated the effect of hub genes on immune cell infiltration in PBC by comparing hub gene expression and immune cell infiltration levels.

Mendelian randomization

All data in this study were obtained from publicly available databases. Two-sample MR was employed to investigate potential causal associations between hub genes and PBC risk, using SNPs as IVs. SNPs associated with PBC were identified using a significance threshold of P < 1e-5, and a linkage disequilibrium cutoff of r² < 0.1 within 10000 bp was applied to each SNP to ensure independence. Data for hub genes were obtained from publicly accessible GWAS data sources. Data for IL10, CCL5, CCL3, and STAT3 were available. SNPs were harmonized to ensure consistent alleles before conducting a two-sample MR analysis. MR analysis was performed using the “TwoSampleMR” package, and the inverse variance weighted (IVW) method was utilized to assess the relationship between the hub gene levels and PBC risk. Additionally, MR-Egger was employed for sensitivity analysis[13].

RESULTS

Screening for DEGs in PBC and HC groups

To conduct our analysis, we retrieved the PBC dataset (GSE79850) from the GEO database and subsequently conducted a rigorous screening process to identify DEGs specific to the PBC dataset. Our detailed investigation revealed a total of 99 DEGs, segmented into 49 genes displaying increased expression and 50 genes showing decreased expression in the PBC group, as compared to the HC group (Figure 1).

Figure 1 Differentially expressed genes between the primary biliary cholangitis and healthy control groups in the GSE79850 dataset. A: A volcano plot was generated to visualize the differential expression analysis results. Genes that were significantly upregulated in primary biliary cholangitis are depicted in red, genes that were significantly down-regulated are depicted in green, and genes that did not show significant differential expression are depicted in black; B: A heatmap was constructed to illustrate differential expression patterns for genes in the GSE79850 dataset. PBC: Primary biliary cholangitis; DEGs: Differentially expressed genes; GO: Gene ontology; KEGG: Kyoto encyclopedia of genes and genomes.

Identification of the PBC-associated module by WGCNA

To gain insight into the potential genetic modules underlying PBC, we utilized WGCNA, employing candidate genes derived from the PBC-related dataset (GSE79850), as illustrated in Figure 2A. Using this method, we successfully categorized the genes into four distinct modules (Figure 2B). Furthermore, we evaluated positive correlation coefficients and identified the "turquoise" module as the most relevant in the context of the GSE89632 dataset (Figure 2C). This specific module shows considerable potential for elucidating its pivotal role in PBC onset and progression.

Figure 2 Identification of primary biliary cholangitis-associated gene modules in the Gene Expression Omnibus dataset using weighted gene co-expression network analysis. A: A dendrogram was generated by clustering genes in the GSE79850 dataset using a topological overlap matrix (1-TOM). Each branch represents a gene, with co-expression modules visualized in different colors; B: A module-trait heatmap displayed the correlation between gene modules and primary biliary cholangitis (PBC) in the GSE79850 dataset. Each module is associated with a correlation coefficient and p-value, indicating the strength and significance of the correlation; C: A scatter plot illustrated the correlation between the turquoise module and PBC in the GSE79850 dataset. The turquoise module exhibits the strongest positive correlation with PBC.

GO/KEGG analyses of key genes

We performed a comprehensive analysis and identified 71 overlapping genes from the WGCNA analysis and DEGs as key genes that may significantly affect PBC development and progression (Figure 3A). We utilized both GO and KEGG frameworks to gain a deeper insight into the functional roles of these 71 overlapping genes. GO enrichment analysis revealed their primary functions in biological processes, encompassing the positive regulation of cytokine production, cell activation, cytokine-mediated signaling pathways, and augmentation of leukocyte activation, as depicted in Figure 3B and C. KEGG analysis highlighted their involvement in complex pathways, such as cytokine-cytokine receptor interactions, Th1 and Th2 cell differentiation, mechanisms underlying inflammatory bowel disease, Th17 cell differentiation pathways, and viral protein interactions with cytokine and cytokine receptor systems, as illustrated in Figure 3D.

Figure 3 Identification and validation of primary biliary cholangitis-associated candidate hub genes in the GSE79850 dataset. A: A Venn diagram was utilized to illustrate the overlap of 115 candidate hub genes; B: GO enrichment analysis of candidate hub genes revealed the functional categories and biological processes associated with primary biliary cholangitis (PBC); C: Kyoto Encyclopedia of Genes and Genomes pathway analysis highlighted the molecular pathways involved in PBC. GO: Gene Ontology; DEG: Differentially expressed gene; WGCNA: Weighted gene co-expression network analysis.

PPI network analysis of overlapped key genes

To delve into the intricate interplay among key genes, we utilized the STRING online platform to generate a protein-PPI network, as depicted in Figure 4A. This network visually captures the complex relationships among these genes, providing insights into their potential functional collaborations in the context of PBC. Additionally, we leveraged Cytoscape software to highlight the top 5 genes that emerged prominently within the hierarchy of the PPI network, providing a focused perspective, as shown in Figure 4B. The top 5 genes with the highest degree values, CD247, IL10, CCL5, CCL3, and STAT3, were identified as hub genes in the network. The color intensity represents degree values, with darker colors indicating higher degree values, thereby highlighting the most significant genes in the network.

Figure 4 Protein-protein interaction network construction. A: A protein-protein interaction (PPI) network was constructed using the identified overlapping candidate hub genes; B: Hub genes in the PPI network were determined using the degree SS algorithm.

Development of a nomogram model for PBC risk prediction

We developed a customized nomogram model, as depicted in Figure 5A, specifically tailored for accurate prediction of PBC risk. This model demonstrates remarkable competence in forecasting the emergence of this condition. Furthermore, ROC curves were generated to assess the diagnostic performance of five crucial genes, namely CD247, IL10, CCL5, CCL3, and STAT3. The AUC values, presented in Figure 5B, highlight the exceptional ability of these genes to distinguish between PBC patients and healthy individuals. The AUC values for CD247, IL10, CCL5, CCL3, and STAT3 were 0.891, 0.969, 0.969, 0.914, and 0.898, respectively.

Figure 5 Prediction of primary biliary cholangitis risk using a nomogram. A: A nomogram model incorporating hub genes was developed to predict the primary biliary cholangitis risk; B: Receiver operating characteristic curve curves were used to evaluate the diagnostic performance of each hub gene. AUC: Area under the receiver operating characteristic curve.

Assessment of immune cell infiltration and the relationship between hub genes and immune cells in PBC

To identify the pertinence of hub genes within the immune microenvironment, the ratio of immune cell infiltration was analyzed using the CIBERSORT algorithm. The findings revealed distinct patterns of immune cell infiltration between PBC and HC (Figure 6). In addition, CD247 expression was significantly correlated with the ratio in T follicular helper cells (Figure 7A). IL10 expression was distinctly associated with resting dendritic cells (Figure 7B). CCL5 expression was significantly associated with the numbers of 3 types of immune cells, including M0 macrophages, naive B cells, and monocytes (Figure 7C). CCL3 expression was distinctly associated with monocytes (Figure 7D). Moreover, STAT3 expression was distinctly associated with 3 types of immune cells, including activated mast cells, M1 macrophages, and resting CD4 memory T cells (Figure 7E). These results provide compelling evidence of a direct association between immune cell infiltration and the hub genes in PBC.

Figure 6 Immune cell infiltration in primary biliary cholangitis. A: Heat map showing the relationship between 22 different types of infiltrating immune cells; B: Heat map displaying variations among 22 types of infiltrating immune cells between primary biliary cholangitis (PBC) specimens and healthy samples; C: Violin plot displaying variations among 22 types of infiltrating immune cells between PBC specimens and healthy samples.

Figure 7 Pertinence between CD247, IL10, CCL5, CCL3, STAT3, and immune cells of primary biliary cholangitis. A: CD247; B: IL10; C: CCL5; D: CCL3; E: STAT3.

There was no causal association between hub genes and the risk of PBC

All SNPs were not considered weak IVs. We examined the causal association between hub genes and PBC using the IVW method. Our findings revealed no significant association between hub genes and PBC risk (Figure 8); IL10 [odds ratio (OR) 1.08, 95%CI: 0.99-1.18, N_SNPs = 6, P = 0.10)], CCL5 (OR 0.95, 95%CI: 0.87-1.04, N_SNPs = 24, P = 0.30), CCL3 (OR 1.00, 95%CI: 0.96-1.05, N_SNPs = 22, P = 0.96), and STAT3 (OR 0.99, 95%CI: 0.91-1.08, N_SNPs = 29, P = 0.82; Figure 8).

Figure 8 Mendelian randomization study results show no causal relationships between the hub genes and primary biliary cholangitis. These include IL10, CCL5, CCL3, and STAT3, as demonstrated using forest plots. A: IL10; B: CCL5; C: CCL3; D: STAT3. MR: Mendelian randomization.

DISCUSSION

Advancements in molecular biology methodologies have sparked heightened interest in exploring genes linked to PBC. By analyzing genome-wide expression patterns, scientists can delve into intricate PBC-associated regulatory networks. Integrating genome-wide expression data enables researchers to identify key genes and dysregulated pathways in PBC, which can then be used as targets to develop therapies that specifically address the underlying molecular disease-related abnormalities. Advancements in molecular biology techniques have provided researchers with powerful tools to study gene roles in PBC progression. The analysis of gene expression profiles and integration of comprehensive genome-wide expression datasets has enabled researchers to decipher intricate regulatory mechanisms underlying PBC, ultimately aiming to identify promising therapeutic avenues for effective disease management.

Our investigation utilized the WGCNA and PPI methodologies to identify pivotal genes related to PBC. This analysis revealed five hub genes-CD247, IL10, CCL5, CCL3, and STAT3-that emerged as highly pertinent to PBC pathogenesis. Subsequently, we constructed a nomogram model integrating multiple factors, including the expression levels of the five hub genes, to accurately predict PBC risk, thereby providing a comprehensive tool for assessing PBC risk. The AUC analysis of the five pivotal genes underscored their proficiency in distinguishing PBC from the control cohort, emphasizing their potential as valuable diagnostic and therapeutic markers. Our nomogram model also emerged as a robust tool for PBC diagnosis and management. Integrating gene expression data and developing predictive models can contribute to more accurate and personalized approaches in PBC diagnosis and treatment.

Using the Cytoscape software, we conducted a PPI network analysis and identified CD247 as the top-ranked gene in the network. CD247, also known as the CD3 zeta chain, is a crucial hub gene associated with PBC. CD247, a transmembrane protein that is an integral component of the CD3 complex, is essential for T-cell receptor signaling and is primarily expressed on T cells. CD247 plays a crucial role in the activation and regulation of T-cell responses. It functions as a signaling subunit within the T cell receptor complex to transmit signals from T cell receptors to intracellular signaling pathways[14,15]. In PBC, CD3 has been implicated in the dysregulation of immune responses and the development of autoimmune liver damage. CD3 is reportedly upregulated in PBC patients, suggesting its involvement in disease pathogenesis[16]. T cells play a significant role in PBC pathogenesis, with abnormal T cell activation and function potentially leading to liver inflammation and damage. In PBC, both CD4+ and CD8+ T cells exhibit quantitative and functional abnormalities. CD4+ T cells in PBC patients reportedly have aberrant regulatory functions, leading to a disrupted immune tolerance and an enhanced autoimmune response[17]. CD8+ T cells are vital in PBC as they infiltrate the liver and have cytotoxic effects on bile duct epithelial cells, resulting in ductal damage[18]. Furthermore, PBC patients have an abundance of natural killer T cells in the liver, which exhibit abnormal activation and function. These natural killer T cells contribute to PBC pathogenesis by producing inflammatory mediators and activating other immune cells[19]. In-depth investigations are imperative to unravel the intricate interplay between T cells and PBC, thereby laying a solid theoretical foundation for the conception of innovative therapeutic approaches.

Interleukin-10 or IL10 plays a pivotal role in modulating immune responses and mitigating inflammation. Originating from diverse immune cells, such as T and B lymphocytes, macrophages, and dendritic cells, IL10 achieves its regulatory functions by engaging with specific receptors on target cells, subsequently triggering a cascade of intracellular signaling events[20,21]. The role of IL10 in PBC has been the subject of extensive research. Some studies have reported decreased levels of IL10 in PBC patients, suggesting that IL10 deficiency could play a role in disease pathogenesis[22,23]. IL10 deficiency may contribute to the dysregulation of immune responses and the development of autoimmune liver damage in PBC. Furthermore, IL10 has been implicated in the modulation of fibrosis and inflammation. IL10 might also have anti-fibrotic and anti-inflammatory effects in the liver, potentially protecting against liver damage and disease progression[24,25]. Despite these findings, the precise mechanisms by which IL10 affects PBC initiation and progression remain elusive. Hence, further research is essential to assess the precise contributions of IL10 in the pathological processes of PBC and to evaluate its therapeutic potential as a target for disease management.

This study marks the first exploration of the potential causal relationship between five hub genes and the risk of developing PBC using a two-sample MR analytical strategy. Utilizing extensive GWAS data for both hub genes (as exposures) and PBC (as the outcome), the MR investigation did not yield evidence indicating a causal relationship between these hub genes and PBC risk. MR methodology has a rigorous design similar to that of prospective randomized controlled trials and serves as a robust tool to counteract the systemic biases commonly encountered in traditional observational research, including confounding factors and the possibility of reverse causality. By incorporating highly precise genotyping methodologies, this approach effectively mitigates the issue of regression dilution due to detection inaccuracies. By leveraging genetic variants as IVs, MR analysis provides a robust framework to assess the causal relationships between exposures and outcomes. However, several limitations are associated with MR studies, including the assumption that genetic variants affect the outcome solely through the exposure of interest and the potential for horizontal pleiotropy. Further research is needed to explore the complex interplay between hub genes and PBC development. The integration of MR analysis in this study contributes to a better understanding of the potential causal relationships between hub genes and PBC risk, providing valuable insights for future investigations and potential therapeutic interventions.

Despite the significant findings, our study has a few other limitations. First, the bioinformatics-driven investigation into hub genes and their potential roles in PBC pathogenesis constitutes an initial step. To substantiate the specific mechanisms of these identified hub genes, subsequent biological experiments, encompassing both in vitro and in vivo studies, are imperative. These experiments would provide additional evidence of the functional roles of hub genes in PBC pathogenesis. Furthermore, it is crucial to note that our MR investigation was solely conducted among individuals of European descent, potentially limiting the broad applicability of our conclusions to diverse populations characterized by differing genetic profiles and environmental exposures. Replicated studies involving diverse populations would be valuable to assess the consistency and robustness of the observed associations. Addressing these limitations in future studies can build upon our findings and contribute to a more comprehensive understanding of the molecular mechanisms underlying PBC development and progression.

CONCLUSION

Using WGCNA-based co-expression network analysis, we identified five hub genes strongly associated with PBC. This discovery provides a strong foundation for future research, enhancing our understanding of the intricate molecular mechanisms underlying PBC. It could propel the development of targeted therapeutic approaches and personalized treatment plans, advancing early diagnosis and optimized PBC management. Ultimately, these advancements are promising for improving patient outcomes and alleviating the overall burden of this liver disease.

ACKNOWLEDGEMENTS

We are grateful to all the participants and investigators involved in the TCGA and GWAS studies included in this study.