3.1MR analysis
Following the MR analysis, the specific results are presented in Table 1. Mannose and TNF-β exhibit negative causality for PBC, whereas Mannose shows positive causality for TNF-β (P < 0.05). Furthermore, neither Mannose nor TNF-β displays horizontal pleiotropy or heterogeneity (P > 0.05).
Table 1
Results of MR analysis of Mannose and TNF-β with PBC
Exposure | Outcome | Method | Nsnp | Beta | Se | Pavle | Pleiotropy | Heterogeneity |
Mannose | PBC | IVW | 10 | -1.621 | 0.699 | 0.020 | 0.830 | 0.113 |
TNF-β | PBC | IVW | 5 | -0.763 | 0.198 | < 0.001 | 0.772 | 0.150 |
Mannose | TNF-β | IVW | 20 | 0.405 | 0.185 | 0.028 | 0.897 | 0.880 |
The results of IVW method analysis showed that the OR was 0.1977 (95% CI: 0.0502 ~ 0.7776,p = 0.020) for PBC as the outcome with Mannose as the exposure factor, 0.4661 (95% CI: 0.3160 ~ 0.6875,p < 0.001) for PBC as the outcome with TNF-β as the exposure factor, and 1.4993 (95% CI: 1.0438 ~ 2.1537,p = 0.028) for Mannose as the exposure factor and TNF-β as the outcome. The remaining methods are shown in Fig. 2.
All IVW estimates in the forest plot were close to zero, indicating the robustness and reliability of the results (Fig. 3A). The distribution of causal effects in the funnel plot was essentially symmetric, with minimal bias, suggesting negligible differences between the instrumental variables (Fig. 3B). The leave-one-out method was employed for sensitivity analysis, excluding each SNP one at a time, recalculating the effect values for the remaining SNPs, and observing the results. This approach tested whether the causal relationship was influenced by a single instrumental variable. The results, compared with MR analysis, consistently showed estimates close to zero, indicating robust MR analysis outcomes and minimal influence on causal correlation by excluding any SNP (Fig. 3C). The slopes of the scatter plots for Mannose and TNF-β in relation to PBC were negative, indicating that both Mannose and TNF-β act as protective factors for PBC. In contrast, the positive slope of the scatter plot for Mannose with TNF-β suggested that Mannose facilitates TNF-β expression (Fig. 3D).
3.2RNA-seq
The quality control results of the scRNA dataset (Fig. 4A) and the downscaling de-batching outcomes (Fig. 4B) were satisfactory, with an inflection point set at 10 (Fig. 4C), categorizing cells into 18 clusters (Fig. 4D). Cell types in the scRNA dataset were automatically annotated as hepatocytes, monocytes, NK cells, B cells, and T cells (Fig. 4E), and the proportions of these cell types in the control and PBC groups were visualized (Fig. 4F). The proportion of hepatic cells decreased sharply, while the proportion of each immune cell increased compared to the control group.
In the scRNA dataset, dot plots and violin plots for LTA clustering were generated, revealing predominant expression in T cells (Fig. 5A, B). Subsequently, the scRNA-TB dataset was divided into MRhigh and MRlow groups based on whether the MR score exceeded the median of all its samples' MR scores. Dot plots and violin plots for LTA clustering were then generated for these groups, revealing predominant clustering in the MRhigh group (Fig. 5C, D).
In the scRNA dataset, GSEA analysis of differential genes in the MRhigh and MRlow groups revealed that the MRhigh group was significantly enriched in pathways related to anti-bacterial and anti-microbial humoral responses, defense against bacteria, plasma lipoprotein particle clearance, and humoral immune responses (Fig. 6A). Meanwhile, the MRlow group showed enrichment in pathways involving αβ T cell activation, proliferation, differentiation, regulation of T cell activation, T cell receptor signaling pathway, and IL-12 synthesis and secretion (Fig. 6A).
Analysis of transcription factor enrichment in the scRNA dataset showed that the MRhigh group was predominantly enriched for BHLHE40 and FOXP1 (Fig. 6B). Meanwhile, the MRlow group exhibited enrichment in transcription factors such as MEF2C, GABPA, TCF7, FOXK2, and MEF2A (Fig. 6B).
The cellular communication results of the scRNA dataset indicated that the TNF signaling pathway was mainly enriched among NK cell-MRlow, Monocyte-MRlow, MRhigh-MRlow, and MRlow-MRlow (Fig. 7A); Monocyte cells were the stronger signal receiver, and the MRhigh and MRlow signals were mostly transduced in NK cells and Monocyte cells (Fig. 7B). Further analysis of the ligand-receptor pairs revealed the presence of NK cell-MRlow, Monocyte-MRlow, MRhigh-MRlow, MRlow-NK cell, MRlow-MRlow, MRhigh-NK cell, MRhigh-MRlow, and MRhigh-MRlow within the TNF-TNFRSF1A and TNF-TNFRSF1B mediated TNF signaling pathway(Fig. 7C).
Moreover, an association analysis of transcription factors downstream of the ligand-receptor pairs indicated that the key transcription factors of the TNF signaling pathway included CREB1, ATF2, ELK4, NFATC1, and RELB for TNF-TNFRSF1A-associated transcription factors (Fig. 8A); the key transcription factors IRF1 and JUN of the TNF signaling pathway among MRlow-MRlow were TNF-TNFRSF1B-associated, and JUN was also a TNF-TNFRSF1A-associated transcription factor(Fig. 8B).
The scRNA-TB group was re-dimensioned and re-clustered, and the results showed that it could be further divided into eight subpopulations (Fig. 9A). Based on the results of a vector-fit temporal analysis, an overview of the progression of the T cell and B cell subpopulations was constructed (Fig. 9C, D), with the pseudo-temporal order of the cells indicated by a color gradient (from red to blue), and the blue rectangles highlighting the starting points in the UMAP cell region. The results indicated that the developmental trajectory progresses from subpopulation 0 to subpopulations 1, 3, and 2. Labeling MR receptor scores alone (Fig. 9B) revealed a decrease in MR receptors over time.
After the bulk dataset was sorted and de-batched (Fig. 10A, B), the differential genes were compared with single-cell MR score differential genes to identify a total of four key genes: DERL3, FKBP2, SSR4, and PRDX4. The genes were screened by applying LASSO, Random Forest, REF-SVM, and GMM machine learning methods to plot their ROC curve (Fig. 10C), and the genes screened by the four machine learning methods were displayed as Wayne diagrams (Fig. 10D) to obtain a core gene, FKBP2. FKBP2 was further validated by applying advanced machine learning methods: Cat Boost, Ng Boost, and XG Boost, and plotted on its ROC curve (Fig. 10E), all showing positive results.