All liver tissue and peripheral blood samples used in this study were collected from the Second Affiliated Hospital of Kunming Medical University, with a collection period from January 2020 to December 2024. The samples included 40 liver biopsy samples from patients with PBC and 40 liver biopsy samples from non-PBC patients. In the PBC patient group, subgroup analysis was further performed according to pathological staging. Additionally, peripheral blood samples were obtained from 15 PBC patients and 15 healthy controls. All samples were approved by the Ethics Committee of the Second Affiliated Hospital of Kunming Medical University (Approval Number: Review-PJ-Department-2023-296), and informed consent was obtained from all patients. All samples were stored and managed in compliance with the hospital’s biobank regulations. The inclusion criteria were as follows: age ≥ 18 years; meeting the EASL 2017 diagnostic criteria for primary biliary cholangitis, including at least two of the following indicators: elevated ALP (≥1.5 times the upper limit of normal), positive AMA or AMA-M2 antibodies, and liver tissue pathology showing small bile duct damage and portal inflammatory infiltration. The exclusion criteria were as follows: coexisting viral hepatitis, autoimmune hepatitis, tumors, other severe systemic diseases, and pregnancy, as well as other pathological conditions unrelated to the study. All samples were obtained from patients at the Second Affiliated Hospital of Kunming Medical University, ensuring the representativeness and consistency of the experimental samples.

Supplementary Table 1 provides details of the 15 PBC patients and 15 normal blood samples used for transcriptome sequencing to create the training set. The GSE119600 dataset (based on the GPL10558 platform) was also downloaded from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). This dataset comprises whole blood samples from 90 PBC patients and 47 controls and was used as the validation set for this study. A total of 1,949 MRGs were downloaded from MitoCarta 3.0 database (8).

Data used for sequencing were evaluated by FastQC package (version 0.11.9) (9). Before analysis, raw sequencing data underwent quality control preprocessing to eliminate adapter sequences and low-quality reads. The resulting clean reads were subsequently aligned to the murine reference genome using HISAT2 (v2.2.1) (10) with standard parameters. Transcript and gene expression quantification were performed using StringTie’s FPKM normalization method. This approach accounted for both mapped read counts and transcript length variations across samples. Expression pattern disparities were subsequently visualized through principal component analysis (PCA). The DESeq2 package (version 1.36.0) (11) was used to identify DEGs between PBC and control samples. Screening criteria were p < 0.05 and |log₂FoldChange (FC)| > 1. Then ggplot2 (version 3.3.5) (12) and pheatmap package (version 1.0.12) (13) were used to draft volcano map and heatmap to display DEGs.

For further accessing genes associated with PBC, WGCNA was conducted using the WGCNA package (version 1.70-3) (14). Initially, hierarchical clustering of all samples was executed to investigate outliers and determine whether their exclusion was necessary to ensure subsequent analyses’ accuracy. Network construction parameters were optimized by selecting a β value that satisfied dual criteria (1): scale-free topology fit (R² ≥ 0.85) and (2) minimal average connectivity. Gene relationship analysis then proceeded through pairwise similarity computation, culminating in phylogenetic tree construction via average linkage hierarchical clustering. Module detection was performed through dynamic tree cutting algorithm implementation, enforcing a minimum module size of 100 genes. Modules demonstrating the strongest association (highest correlation coefficients) with phenotypic traits were selected as pivotal modules, whose member genes were classified as biologically relevant candidates.

By taking the intersection of DEGs, key module genes and MRGs, candidate genes were obtained. To find potential functions and pathways enriched by candidate genes, GO and KEGG enrichment analyses were performed by clusterProfiler package (p < 0.05 and count ≥ 1) (version 4.0.2) (15).

LASSO and SVM-RFE were carried out to screen signature genes. Interaction of signature genes derived from 2 machine-learning algorithms was identified as candidate hub genes. Receiver operating characteristic (ROC) curves were then plotted for all samples in the training and validation cohorts using the pROC package (version 4.0.2) to evaluate the ability of the candidate hub genes to distinguish between the PBC patient group and the control group (16). Genes with an area under the curve (AUC) of over 0.7 in both the training and validation sets were selected as the final hub genes and their expression levels were analyzed. Furthermore, correlation between hub genes was explored by spearman correlation analysis. Cell-PLoc 2.0 performed subcellular localization on proteins of hub genes (17).

In order to explore the biological function of hub genes, GSEA was carried out. Background gene set c5.go.v7.4.entrez.gmt and c2.cp.kegg.v7.4.entrez.gmt were downloaded from GSEA website. Samples were divided into high and low expression groups based on expression median value of hub genes. GSEA of all genes in 2 expression groups was processed with threshold |NES| > 1, NOM P < 0.05 and q < 0.25. Back-propagation neural network (BPNN) was constructed by neuralnet package (version 4.0.2) (18) to evaluate accuracy of hub genes.

The similarity of GO terms of hub genes was calculated by GOSemSim package (version 2.18.1) (19).

To analyze differences in the immune microenvironment between the PBC and control groups, CIBERSORT (version 0.1.0) was used to evaluate the infiltration levels of 22 immune cell subpopulations (20) Gene and immune cell expression matrices were processed using a deconvolution algorithm to obtain an immune infiltration profile, visualized using ggplot2 (version 3.3.5). To identify immune cells showing significant differences between the two groups, we analyzed differences in immune cell composition between samples using the Wilcoxon test, setting p < 0.05 as the significance threshold. The results were visualized using ggplot2 (version 3.3.5).

MicroRNAs (miRNAs)-targeting hub genes predicted through miRDB (http://mirdb.org) and miRWalk database. Key miRNAs were acquired by overlapping miRNAs obtained from the 2 databases (energy < -30). Then, lncRNAs-targeting miRNAs were obtained from the miRTarBase database. Finally, Cytoscape software was used to visualize the lncRNAs-miRNA-mRNA network. To further understand the regulatory effects of hub genes in disease, TFs-targeting hub genes were predicted using Cistrome database (RF score > 0.75), and the TF-mRNA-miRNA network was constructed to visualize potential relationships.

Correlated diseases of hub genes were analyzed by DisGeNET database (http://www.disgenet.org/), and the hub genes-diseases co-expression network was visualized by NetworkAnalyst (http://www.networkanalyst.ca). Based on hub genes, small molecule drugs were found by DGIdb.

From January 2018 to October 2024, a total of 30 paraffin-embedded liver tissue samples were collected from patients diagnosed with primary biliary cholangitis (PBC) at the Department of Pathology, The Second Affiliated Hospital of Kunming Medical University. All diagnoses were confirmed by pathological examination. For comparison, 17 normal liver tissue samples were used as controls. PBC patients were stratified into four histological stages according to Ludwig’s classification. Immunohistochemical staining was performed to assess the protein expression of SHANK2 and TGM2. The staining protocol strictly followed the instructions provided by the antibody manufacturers. Sections were deparaffinized, rehydrated, and subjected to antigen retrieval, followed by incubation with primary antibodies against SHANK2 and TGM2, respectively. Signal detection was achieved using an HRP-conjugated secondary antibody and DAB chromogen. Expression levels were semi-quantitatively assessed based on staining intensity and percentage of positive cells. Comparisons between groups and different stages were conducted using appropriate statistical methods, with p-values < 0.05 considered significant.

Peripheral blood samples were obtained from 10 patients with PBC and 8 healthy controls. Peripheral blood mononuclear cells (PBMCs) were isolated using Ficoll density gradient centrifugation. Total RNA was extracted from PBMCs using TRIzol reagent, and the RNA concentration and purity were determined using a Nanodrop spectrophotometer. Complementary DNA (cDNA) was synthesized using a reverse transcription kit, and qPCR was performed using SYBR Green Master Mix on a real-time PCR system. The relative mRNA expression levels of SHANK2 and TGM2 were quantified using the 2^−ΔΔCt method and normalized to GAPDH as an internal control. Statistical comparisons between groups were conducted using a t-test, and a p-value of < 0.05 was considered statistically significant.

All statistical analyses were conducted in R software and Microsoft Excel. The Wilcoxon test and t test were utilized to assess the differences between different groups. The P value less than 0.05 was considered statistically significant.

Alignment rate of 30 sequencing samples was all above 90% through alignment analysis and base error rate was less than 1/1000 ( Supplementary Table 2 ). And FPKM distribution of each samples showed that sample sequencing effect was decent ( Figure 1A ). PCA diagram and Boxplot showed that overall sequencing quality was high ( Figures 1B-C ). Furthermore, 1,905 DEGs were identified between PBC and control samples, of which 1,706 up-regulated and 199 down-regulated ( Figures 1D-E ).

We performed WGCNA to identify key module genes associated with PBC. Based on hierarchical clustering of all sequencing samples, one outlier (YQKM2B14) was detected and excluded from downstream analyses ( Figure 2A ). Subsequently, we selected the optimal soft-thresholding power (β = 9), which yielded a scale-free topology fit index R² greater than 0.85 with mean connectivity close to zero, indicating that the constructed network closely approximated a scale free topology ( Figure 2B ). Using this parameter, a gene co-expression network was established and a total of 10 distinct modules were identified, each containing at least 100 genes ( Figure 2C ). Lastly, Module–trait relationship analysis revealed that the MEroyalblue module showed the strongest positive correlation with PBC(|cor| = 0.44, P < 0.05), and 257 genes within this module were considered key module genes ( Figures 2D, E ). Byoverlapping the 1,905 DEGs, the 257 MEroyalblue module genes, and 1,949 mitochondria-related genes (MRGs), weidentified seven candidate genes: AURKA, C1QC, MAOA, RAD51, SHANK2, SPARCL1, and TGM2 ( Figure 2F ; Supplementary Table 3 ). To further explore the biological roles of these candidate genes, GO and KEGG enrichment analyses were performed. In total, 4 biological process (BP), 12 cellular component (CC), 17 molecular function (MF), and 13 KEGG pathways were significantly enriched. The top representative terms included histone modification, synapse organization, oxidoreductase activity, and tyrosine metabolism Figure 2G , These findings indicate that the candidate genes are functionally linked to epigenetic regulation, energy metabolism, and amino acid metabolism, suggesting their potential involvement in the pathogenesis of PBC.

LASSO regression algorithm was used to select MAOA, SHANK2 and TGM2 as signature genes (lambda.min = 0.08) ( Figure 3A ). Then SVM-RFE algorithm screened 6 signature genes2 (TGM2, AURKA, C1QC, MAOA, SHANK2 and RAD51) when maximal accuracy was 0.633, minimal RMSE was 0.367 ( Figure 3B ). Intersection of signature genes (MAOA, SHANK2 and TGM2) was taken and identified as candidate hub genes ( Figure 3C ). To evaluate these hub genes’ ability to distinguish PBC from control samples, ROC curves were plotted for all samples in the training and validation sets. The results showed that the AUC for both SHANK2 and TGM2 exceeded 0.7, whereas the AUC value for MAOA was 0.607 in the training set. Thus, SHANK2 and TGM2 were obtained as hub genes for subsequent analysis ( Figures 3D, E ). Expression of SHANK2 and TGM2 was higher in PBC group ( Figure 3F ).

Spearman correlation analysis was used to determine whether there was correlation between hub genes. Results indicated that SHANK2 and TGM2 had a significant positive correlation (r = 0.42, p < 0.05) ( Figure 4A ). Similarity analysis showed that the similarity score of SHANK2 and TGM2 was greater than 0.4, which indicated that hub genes had strong similarity ( Figure 4B ). Subcellular localization showed TGM2 was mainly expressed in the plasma membrane, while SHANK2 was mainly expressed in nucleus and plasma membrane ( Figure 4C ). Based on hub genes, a BPNN was constructed and prediction correlation and error respectively were 0.8199 and 0.1326, indicating a small prediction bias of the BPNN model ( Figure 4D ). Through GSEA, SHANK2 correlated with 921 GO terms and 56 KEGG pathways which mainly were B cell receptor signaling pathway, T cell receptor signaling pathway and wnt signaling pathway ( Figure 4E ). TGM2 enriched in 1364 GO terms and 74 KEGG pathways, including oxidative phosphorylation, pyrimidine metabolism and purine metabolism ( Figure 4F ).

The percentage abundance of immune infiltration cells for PBC and control samples is shown in Figure 5A . Immune profiling identified five differentially abundant leukocyte subsets between cohorts. The control group demonstrated CD8+, dormant NK cells, and M2 macrophages, contrasting with PBC specimens, which displayed enrichment of central memory CD4+ T cells and undifferentiated M0 macrophages ( Figure 5B ). In these immune cells, activated NK cells positively correlated with activated mast cells (r > 0.4, p < 0.05), while monocytes were negatively related with neutrophils (r < -0.4, p < 0.05) ( Figure 5C ). Finally, hub genes correlated with differential immune cells. SHANK2 was positively related to M0 and M1 macrophages, while TGM2 had significant positive correlation with monocytes ( Figure 5D ).

In order to explore involved regulatory mechanism of hub genes, lncRNA-miRNA-mRNA network was constructed with 45 key miRNAs and 21 lncRNAs. SLFN12L, COLCA1 and LINC00598 regulated SHANK2 through hsa-miR-665, and DNAH10OS, SLFN12L, TSPEAR-AS2 and HTR5A-AS1 regulated TGM2 by hsa-miR-6720-5p ( Figure 6A ). Then 17 TFs were predicted, and TF-mRNA-miRNA network was constructed. In TF-mRNA-miRNA network, TFAP2C, REST, PPARG, etc regulated TGM2, while only TP53 regulated SHANK2 (Figure 6B ). Based on hub genes, 49 diseases were predicted by DisGeNET, bipolar disorder was disease predicting by SHANK2 and TGM2 simultaneously ( Figure 6C ). In addition, 41 small molecule drug types were predicted by hub genes, 39 types-targeting TGM2 (CISPLATIN, VITAMIN A, etc) and 2 types-targeting SHANK2 (calcium, N-Methyl-D-aspartic acid (NMDA)) were included ( Figure 6D ).

To further validate the bioinformatics findings at the protein level, immunohistochemistry was employed to evaluate SHANK2 and TGM2 expression in liver tissues. SHANK2 was significantly overexpressed in the liver tissues of PBC patients compared to the control group, with prominent expression differences observed in stages I, II, and III of PBC (P < 0.05) ( Figure 7A ). Similarly, TGM2 protein expression was markedly elevated in the PBC group, with statistically significant increases noted across all four pathological stages (P < 0.05) ( Figure 7B ). These results suggest that both SHANK2 and TGM2 proteins are upregulated in the progression of PBC and may be involved in stage-specific pathological processes.

To investigate SHANK2 and TGM2 expression at the transcriptional level, qPCR analysis was performed on RNA extracted from PBMCs of PBC patients and healthy controls. The results demonstrated significantly elevated mRNA expression levels of both SHANK2 and TGM2 in the PBC group compared to controls (P < 0.05) ( Figure 8A, B ). These findings confirm the upregulation of SHANK2 and TGM2 in PBC not only in liver tissue but also in peripheral immune cells, indicating a potential systemic involvement in disease pathogenesis.