Single-cell and bulk transcriptomics with machine learning decode B cell hub genes and diagnostic biomarkers in anal fistula

Abstract

BACKGROUND

Anal fistula (AF) is a perianal inflammatory disorder with a complex etiology and high recurrence rates, causing substantial functional impairment. The immune microenvironment of AF is heterogeneous, but its potential mechanisms of action remain unclear.

AIM

To establish a multimodal framework for AF by integrating single-cell and bulk RNA sequencing datasets.

METHODS

High-dimensional weighted gene co-expression network analysis combined with five machine learning (ML) algorithms were used to identify B-cell-associated hub genes. Immune infiltration profiling, functional enrichment, pseudotime trajectory, and intercellular communication analyses were performed. A gene-based predictive model was constructed and validated using receiver operating characteristic curves and a nomogram to evaluate diagnostic performance.

RESULTS

About 29 cellular clusters with enriched B cells were identified in the AF tissues. Among 129 B-cell-related candidate genes, core genes such as long intragenic noncoding RNA p53-induced transcript, lipopolysaccharide-responsive and beige-like anchor, and spleen tyrosine kinase, were consistently selected for all ML algorithms. These genes showed strong correlations with multiple immune cell types and inflammatory pathways. B cells associated hub genes enrichment was detected in the pathway related to transforming growth factor beta/mitogen-activated protein kinase signaling, primary immunodeficiency, and natural killer cell-mediated cytotoxicity, and these B cells acted as central signaling nodes promoting inflammatory amplification, angiogenesis, and fibrosis. The predictive model demonstrated robust diagnostic accuracy, with an area under the curve of 0.858.

CONCLUSION

This study is the first to integrate single-cell and bulk transcriptomics with ML to systematically decode the B-cell associated immune network in AF. The identified hub genes may serve as crucial diagnostic biomarkers and therapeutic targets.

Key Words: Anal fistula; Single-cell RNA sequencing; Bulk RNA sequencing; Machine learning; Immune microenvironment

Core Tip: This study is the first to combine single-cell RNA sequencing (RNA-seq) with bulk RNA-seq and machine learning to systematically characterize B-cell-associated immune networks in anal fistula (AF). Three hub genes, LINC-PINT, LRBA, and SYK, were identified as central signaling nodes linking inflammatory amplification, angiogenesis, and fibrosis within the AF microenvironment. A predictive model incorporating these genes achieved robust diagnostic accuracy, offering promising biomarker candidates for precision diagnosis and targeted therapeutic intervention in AF.

INTRODUCTION

Anal fistula (AF) is a common perianal disorder with an incidence rate of 12.3 per 100000 in men and 5.6 per 100000 in women[1], suggesting that men are affected approximately twice as frequently as women. This condition typically occurs around 40 years of age. Although the exact etiology and pathogenic mechanisms remain unclear, AF is generally caused by an infection of the anal glands (located within the intersphincteric space). The infection diffuses along the pathway of least resistance forming a persistent tract and an aberrant connection between the anorectal canal and the perianal skin[2]. Clinically, patients often experience recurrent abscesses or poorly draining fistulas that lead to chronic or intermittent symptoms such as pain, purulent discharge, and social discomfort. Managing complex fistulas poses a significant surgical challenge. Despite ongoing advancements in surgical techniques, failure rates remain high, and postoperative anorectal dysfunction and recurrence are prevalent. Relapse occurs in 20%-30% of patients with complex fistulas[3], who may also be prone to an increased risk of sphincter impairment, local fibrosis and scarring, and defecatory disorders, thereby lowering their quality of life[4]. Therefore, accurate early prediction and comprehensive assessment of fistula onset and progression are crucial for developing effective preventive and therapeutic strategies.

Immune dysregulation is a pivotal factor in AF disease pathogenesis and progression. A comprehensive analysis presented by the 5^th Scientific Workshop of the European Crohn’s and Colitis Organization demonstrated that in AF, a marked immune cell infiltration was observed including T cells, macrophages, CD20⁺ B in the local inflammatory microenvironment[5]. This highlights the role of B cells in chronic inflammation and tissue remodeling characteristic to AF. Protein-protein interaction network studies have further identified core genes shared between AF and colorectal cancer[6] including CDKN2A and TIMP1[7], suggesting that dysregulation of immune-associated gene expression may play an important role in disease development. Aberrant expression of CDKN2A is associated with disrupted cell cycle regulation[8] indicating its potential involvement in AF pathogenesis by disrupting the balance between cellular proliferation and apoptosis. TIMP1 expression is significantly associated with the infiltration levels of B cells, CD8⁺ T cells, and macrophages[9]. During chronic inflammation in AF, diminished TIMP1 expression may impair the negative regulation of these immune cells causing prolonged inflammation, delaying tissue repair, and exacerbating local immune dysregulation and fistula chronicity[10,11].

A previous systematic review and meta-analysis validated the safety and effectiveness of mesenchymal stem cell transplantation for the treatment of complex peri AFs[12]. Notably, mesenchymal stem cells can trigger interleukin-10 expressing regulatory B cells, which play a key role in anti-inflammatory reactions and immune tolerance, effectively ameliorating autoimmune diseases such as experimental colitis and systemic lupus erythematosus[13,14]. Furthermore, dysbiosis of the intestinal microbiota can weaken the mucosal barrier by inducing epithelial-mesenchymal transitions and inflammatory responses[15-17]. As primary mediators of mucosal immunity, B cells can activate, differentiate, and produce antibodies in response to microbial alterations and inflammatory stimuli, amplifying or modulating the local inflammatory environment[18-20]. It is possible that B cells underlie intestinal microbiota population, barrier disruption, and immune responses observed in the pathogenesis of AF. B cells might also contribute to chronic inflammation and impaired fistula healing. A recent single-cell sequencing study demonstrated that B cells undergo class-switch recombination and somatic hypermutation under chronic inflammatory conditions, thereby facilitating plasma cell differentiation and local antibody secretion[21]. Although these processes may enhance initial defense against pathogens, these can also cause tissue damage and fibrosis[22,23], which exacerbate fistula formation and impair wound healing. Therefore, elucidating the expression profiles and functional roles of immune imbalances particularly those involving B cells and their regulatory genes, may provide new mechanistic insights and identify therapeutic targets for the treatment of AF.

A significant strength of this study was the prospective establishment of a rare clinical biobank comprising of an integrated multimodal transcriptomic dataset for AF. This dataset synergistically combines high-resolution single-cell RNA sequencing (scRNA-seq) with bulk RNA sequencing (RNA-seq), both executed on a unified experimental platform (10 × Genomics Chromium v3.0 + DNBSEQ-T7). Building on this foundation, we developed and implemented a high-dimensional weighted gene co-expression network analysis (hdWGCNA) that was innovatively integrated with ensemble machine learning (ML) algorithms. This systematic approach enabled the successful deconvolution of key immune hub gene networks and their regulatory interactions within the complex immune microenvironment (IME) of AFs. Our study identified a novel repertoire of highly translatable biomarker targets, thereby establishing a molecular foundation for advancing precise diagnosis, targeted therapy, and prognostic assessment of AF.

MATERIALS AND METHODS

Study participants and tissue collection

Tissue samples for scRNA-seq analysis were collected from patients diagnosed with cryptoglandular AF (n = 3) and from healthy individuals (n = 3) at the Department of Proctology of China-Japan Friendship Hospital. For bulk RNA-seq analysis (n = 38), fresh tissue specimens were collected intraoperatively including 18 patients with cryptoglandular AF and 20 healthy controls. Diagnosis of cryptoglandular AF was established by proctology research team through a comprehensive evaluation of clinical manifestations, physical examinations, ultrasound imaging, and histopathological criteria. Transrectal ultrasonography revealed hypoechoic or anechoic regions within each lesion occasionally interspersed with hyperechoic signals indicative of gas. Tracing these hypoechoic tracts identified the internal opening which presented as a mucosal defect characterized by focal interruption of mucosal continuity and localized mucosal elevation or depression. Outward tracing showed hypoechoic tracts extending toward the external cutaneous opening. The fistula tracts appeared as single or multiple branches with longitudinal views demonstrating cord-like hypoechoic channels and transverse views displaying round cystic structures; secondary branches were occasionally observed around the primary tract. The AF specimens were reevaluated by skilled pathologists, and the diagnosis was confirmed using hematoxylin and eosin-stained sections (Figure 1).

Figure 1 Histopathology of the internal-opening in cryptoglandular anal fistula (hematoxylin and eosin). A-D: Low-power views of internal-opening tissue; E-H: Higher-power images of the boxed regions. Red dashed contours delineate gland-like structures. A 1-2 cell thick layer of spindle-shaped myoepithelial cells with goblet cells surrounding an anal gland, surrounded by florid chronic inflammatory infiltrates with fibroblast proliferation and collagen deposition. These features are consistent with a glandular origin of the fistula and were reviewed by a senior pathologist.

Tissue dissociation and single-cell suspension preparation

Samples were obtained from patients and rapidly rinsed twice with pre-cooled RPMI 1640 medium containing 0.04% bovine serum albumin (BSA) under sterile conditions. These tissues were then mechanically dissociated into roughly 0.5-mm³ fragments using surgical scissors and transferred into a newly prepared enzymatic digestion solution. The mixture contained RPMI 1640 medium (Corning, Cat. No. 10-040-CVR), 0.04% BSA (MACS, Miltenyi Biotec, Bergisch Gladbach, Germany, Cat. No. 1000076), and 0.2% collagenase II (Gibco, Waltham, MA, United States, Cat. No. 17101015). These samples were cultivated for 30-60 minutes with gentle inversion every 5-10 minutes at 37 °C. The digested cell suspension was screened up to two times using a BD Falcon 40 μm cell strainer (BD, Franklin Lakes, NJ, Cat. No. 352340). The filtrate was centrifuged at 300 × g for 5 minutes at 4 °C. The cell pellet was re-suspended in appropriate medium, blended with an equal volume of red blood cell lysis buffer (Miltenyi Biotec, Cat. No. 130-094-183), and incubated for 10 minutes at 4 °C. After centrifugation at 300 × g for 5 minutes, the supernatant was discarded. The pellet was rinsed once with medium and centrifuged at 300 × g for 5 minutes. The supernatant was discarded. Finally, the cells were re-suspended in 1 mL RPMI 1640 medium (Corning, Cat. No. 10-040-CVR) containing 0.04% BSA. The single-cell suspension concentration and cell activity were evaluated using a Luna-FL cell counter (Logos Biosystems, Anyang, South Korea). The details for RNA isolation and library preparation for bulk RNA-seq are provided as Supplementary material.

RNA isolation and library preparation

Fresh target tissues of suitable dimensions were excised and rinsed with ice-cold enzyme-free water or phosphate-buffered saline to remove blood and debris. Specimens were immediately immersed in liquid nitrogen for rapid freezing, then transferred into pre-cooled cryovials following complete freezing, and subsequently stored in liquid nitrogen prior to long-term preservation at -80 °C. Total RNA was extracted using TRIzol reagent (Invitrogen, Waltham, MA, United States). RNA purity and quantification were performed using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, United States). RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States), and libraries were prepared using the VAHTS Universal V6 RNA-seq Library Prep Kit. Transcriptome sequencing was performed by OE Biotech Co. Ltd. (Shanghai, China).

ScRNA-seq analysis

We used the Seurat package[24] to create a Seurat object from a single-cell expression matrix. The proportion of mitochondrial gene expression was indicative of cellular homeostasis because increased mitochondrial content frequently signals cellular stress. Consequently, cells with mitochondrial gene content exceeding 30%, total counts (nCount) above 30000 or below 200, or detected features (nFeature) > 5000 or < 200, were excluded from the analysis. Sequencing depth was subsequently normalized using the “NormalizeData” function and 2000 highly variable genes were chosen for in-depth analysis. Principal component (PC) analysis[25] was conducted to identify significant components and ElbowPlot was used to visualize their distribution. Fifteen PCs were selected for dimensionality reduction using uniform manifold approximation and projection (UMAP). Ultimately, cells were clustered into 29 subpopulations using the “FindClusters” function, and the resolution parameter was initialized to 0.5. Cell type labeling was performed using canonical marker genes. Specifically, the markers CD3D, CD3E, and NKG7 were indicative of T/natural killer (NK) cells; BANK1 and CD79A were associated with B cells; MZB1 and IGHG1 were characteristic of plasma cells; FCGR3B and CXCL8 identified neutrophils; LYZ and MS4A6A were linked to monocyte-like cells; GATA2 and KIT were markers for mast cells; PECAM1 and VWF were related to endothelial cells; APOD and LUM corresponded to fibroblasts; KRT6A and CDH1 were indicative of epithelial cells; MYLK and ACTA2 were associated with smooth muscle cells; and S100B and CDH2 were markers for glial cells. Regarding batch effects, all samples were processed on a uniform experimental platform (10 × Genomics Chromium v3.0) with consistent library preparation and sequencing protocols (DNBSEQ-T7), effectively controlling technical batch variation at the experimental level. At the analytical level, the Harmony algorithm was applied to integrate the dimensionality-reduced data and eliminate potential technical bias across samples.

Use of high-dimensional WGCNA to identify key B-cell-associated genes

HdWGCNA effectively manages high-throughput single-cell transcriptomic datasets, facilitating the identification of gene expression patterns across diverse cell types. Using hdWGCNA, cell type-specific co-expression networks have been constructed to identify gene modules and co-express genes within these networks[26]. This approach offers a novel perspective for elucidating molecular mechanisms underlying complex diseases. To identify pivotal genes associated with B cells in AF and control groups, we used hdWGCNA via the hdWGCNA R package, employing a soft-thresholding power of β = 5. Additionally, we used the FindMarkers function to determine differentially expressed genes of B cells in the AF and control groups, implementing a minimal expression threshold (min.pct) of 0.25 to ensure robust filtering.

ML algorithms for identifying potential hub genes

In this study, optimal disease-associated feature genes were identified using five ML algorithms: Extreme gradient boosting (XGBoost)[27], random forest (RF)[28], least absolute shrinkage and selection operator (LASSO)[29], Boruta[30], and support vector machine (SVM)[31]. By integrating these algorithms and leveraging their unique feature selection and validation methodologies, we determined optimal feature genes with significant diagnostic or predictive potential.

Gene Ontology/Kyoto Encyclopedia of Genes and Genomes enrichment analysis

Gene Ontology (GO) enrichment analysis[32], a common technique for extensive functional annotation, was employed to assess biological processes, molecular functions, and cell components linked to hub genes. Simultaneously, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis[33] was implemented to examine the signaling and disease-associated pathways. These analyses were conducted using the clusterProfiler R package[34] enabling comprehensive insight into the biological effects of the identified genes.

Gene set enrichment analysis and immune-related analysis

Gene set enrichment analysis (GSEA) is a computational approach used to determine whether a predefined set of genes shows statistically significant concordant variations between two biological states. This algorithm evaluates the distribution patterns of gene sets ranked by their correlation with a specific phenotype, thereby assessing the overall contribution to that phenotype[35]. To explore differences in biological processes between groups, we obtained the “c2.cp.kegg.v7.0” gene set from the molecular signatures database[36]. Using bulk RNA-seq data, we conducted GSEA to assess enrichment of gene sets associated with potential hub genes. Additionally, we implemented single-sample GSEA (ssGSEA) using the GSVA R package[37] to measure the composition and relative abundance of 28 immune cell types in AF samples. We evaluated the associations between the expression levels of potential hub genes and various immune cell populations. Furthermore, we searched for immune-associated functional gene sets in the ImmPort database to examine the links between these key genes and multiple immune signaling pathways thereby providing deeper insights into their immunological relevance.

Development of a gene-based predictive model for AF and visualization using a nomogram

Based on RNA-seq data from patients with AF, receiver operating characteristic (ROC) analysis was performed to individually evaluate the diagnostic performance of key predictive hub genes. For each target gene, the area under the ROC curves (AUC) and statistical significance (P values < 0.05) were calculated. The diagnostic potential of each gene was quantified by the AUC value, with those above 0.7 indicating strong potential as diagnostic biomarkers for AF. To fully assess the combined predictive capacity of these hub genes and enhance individualized risk assessment, we created a multivariable logistic regression model to estimate the probability of disease occurrence. Subsequently, we created a nomogram (also known as a risk-prediction chart) to visually depict individualized risk predictions and facilitate clinical interpretation.

Pseudotime trajectory analysis and cell communication

To clarify the temporal dynamics of pivotal genes in B cells, we conducted a pseudotime trajectory analysis using monocle R package[38]. The gene-cell expression matrix extracted from the Seurat-derived B-cell subset served as input, with default parameters employed to reconstruct developmental trajectories. To comprehensively explore the effects of these genes in B cells from patients with AF, we used the CellChat R package to quantitatively infer and analyze intercellular interaction networks based on scRNA-seq data[39]. Interaction landscapes between cell populations were observed using circle plots whereas bubble plots were generated to highlight critical ligand-receptor pairs mediating intercellular signaling.

Isolation of B cells and reverse transcription-quantitative PCR analysis

B cells were isolated from fresh tissue samples using fluorescence-activated cell sorting. Briefly, single-cell suspensions were labeled with fluorochrome-conjugated antibodies targeting CD3 and CD19. To ensure high purity and eliminate T-cell contamination, B cells were precisely identified and sorted using a CD3^-CD19⁺ gating strategy. After sorting, total RNA was extracted from the purified B cells and the mRNA expression levels of the target genes were quantified by reverse transcription-quantitative PCR (RT-qPCR). Comprehensive protocols detailing tissue dissociation, flow cytometry gating parameters, and primer sequences are provided in Supplementary material.

Statistical analysis

Data processing and statistical analyses were performed using R software (v 4.2.2). For statistical comparisons between the two groups, continuous variables were assessed using the Mann-Whitney U test, also known as the Wilcoxon rank-sum test. All statistical tests were two-tailed, with P values < 0.05 were considered statistically significant.

RESULTS

Histopathological results

Histopathological analysis of the internal opening identified distinctive anal gland structures characterized by the presence of spindle-shaped myoepithelial cells and goblet cells encircled by chronic inflammatory infiltrates. These morphological characteristics, corroborated by a senior pathologist confirmed that the collected specimens were anal gland tissues of glandular origin (Figure 1).

Single-cell analysis of AF

The detailed clinical characteristics of all individuals included in the scRNA-seq analysis are provided in Supplementary Table 1. We analyzed scRNA-seq data derived from patients with AF and controls. Using the UMAP dimensionality reduction, 29 distinct cellular clusters were identified (Figure 2A). The distribution of cells across different groups and individual samples was visualized using UMAP plots (Figure 2B and C). We labeled each cluster according to canonical cell type markers resulting in the identification of 11 major cell types: B cells, neutrophils, glial cells, smooth muscle cells, monocytes, plasma cells, mast cells, T/NK cells, epithelial cells, fibroblasts, and endothelial cells (Figure 2D). We quantified the number of cells within each of these 11 subpopulations and found that fibroblasts were the most abundant (n = 21690), and glial cells were the least represented (n = 361; Figure 2E). We also determined the relative percentages of these cell types in the dataset (Figure 2F). We performed further UMAP clustering analysis, specifically on the B-cell population (Supplementary Figure 1). We identified three major B-cell subsets: Naïve, memory, and germinal center.

Figure 2 Single-cell transcriptomic landscape of anal fistula and control samples. A: Uniform manifold approximation and projection (UMAP) plot showing the clustering of all single cells into 29 distinct Seurat clusters; B: UMAP visualization colored by samples, including anal fistula (AF) and control samples; C: UMAP plot of annotated major cell types across AF and control groups; D: Dot plot displaying the expression of canonical marker genes used to define 11 major cell types. The size of the dots represents the percentage of cells expressing each gene, and the color indicates the average expression level; E: Bar plot showing the absolute number of cells for each of the 11 major cell types; F: Proportional distribution of the 11 major cell types across all six samples. UMAP: Uniform manifold approximation and projection; AF: Anal fistula.

HdWGCNA identified potential hub genes associated with B cells

We created a line chart to illustrate the variations in cell type composition between the AF and control groups. Among the analyzed cell types, B cells demonstrated the most substantial differences with a notably higher proportion in the AF group (Figure 3A). Subsequently, we conducted hdWGCNA to identify the hub genes associated with B cells using a soft-threshold power of 5 (Figure 3B). Fourteen gene modules were identified (Figure 3C and D). Further expression profiling complemented by bubble plot visualization indicated that genes from the green, black, pink, magenta, yellow, green-yellow, and salmon modules were predominantly expressed in B cells (Figure 3E). Differential expression analysis of B cells from the AF and control groups revealed 764 genes that were upregulated in the AF group (Figure 3F) supporting the hypothesis that B cells play an irreplaceable role in the etiology of AF. These findings provide a convincing rationale for the in-depth investigation of B-cell-related molecular alterations.

Figure 3 B cell-specific differences and co-expression modules in anal fistula vs control. A: Line plot showing the distribution of 11 major cell types across anal fistula (AF) and control samples. B cells exhibited the largest difference and were markedly elevated in the AF group; B: Determination of soft-thresholding power in high-dimensional weighted gene co-expression network analysis (hdWGCNA); C: Dendrogram of gene modules identified by hdWGCNA, revealing 14 co-expression modules in B cells; D: Uniform manifold approximation and projection visualization of module gene expression patterns across cells, colored by different module identities; E: Dot plot showing the expression levels of genes from each module across different cell types. Green, black, pink, magenta, yellow, green-yellow, and salmon modules were enriched in B cells; F: Volcano plot showing differentially expressed genes in B cells between AF and control groups. A total of 764 genes were significantly upregulated in the AF group. AF: Anal fistula; hdWGCNA: High-dimensional weighted gene co-expression network analysis.

Functional enrichment analysis

We intersected the differentially expressed genes with those identified from the hdWGCNA modules, and identified 129 overlapping genes that were critical for B-cell function within AF (Figure 4A). KEGG pathway enrichment analysis showed that these genes were significantly enriched in several pathways including actin cytoskeleton regulation, Fc gamma R-mediated phagocytosis, and the BCR signaling pathway (Figure 4B). As suggested by GO biological processes enrichment analysis, these genes were involved in biological processes such as protein-containing complex assembly modulation, nucleobase-containing compound catabolism, and B-cell differentiation (Figure 4C). GO molecular functions analysis indicated enrichment of GTPase regulator activity, nucleoside-triphosphatase regulator activity, and protein serine/threonine kinase activity (Figure 4D). GO cell components analysis revealed significant associations with actin-based cell projections (Figure 4E).

Figure 4 Functional enrichment analysis of candidate genes in anal fistula-related B cells. A: Venn diagram showing the intersection of genes identified from hdWGCNA and differentially expressed genes in B cells, resulting in 129 overlapping genes; B: Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of the 129 overlapping genes; C: Gene Ontology (GO)-biological process enrichment analysis; D: GO-molecular function enrichment analyses; E: GO-cellular component enrichment analysis. hdWGCNA: High-dimensional weighted gene co-expression network analysis.

Identification of hub genes by ML

We conducted a comprehensive analysis of 129 previously identified hub genes using RNA-seq data from AF samples to screen for 19 differentially expressed genes. These genes included STK17B, LINC-PINT, TERF2IP, CYTIP, CMTM6, BLNK, PARP8, DMTF1, VAV3, LRBA, PTPRC, JAM3, EML4, SAMSN1, CYRIB, SYK, CDC42SE2, COP1, and MARCHF3, all of which were upregulated in the AF group (Figure 5A). We evaluated these 19 genes using various ML algorithms to identify predictive features of AF. SVM analysis showed that the model comprising of 13 variables achieved the highest precision and lowest error rate (Figure 5B). LASSO regression highlighted three critical genes, LINC-PINT, LRBA, and SYK (Figure 5C and D). RF analysis provided a ranking and visualization of the 15 most significant genes (Figure 5E). The Boruta analysis identified 6 key genes: BLNK, LINC-PINT, LRBA, MARCHF3, SYK, and TERF2IP (Figure 5F). Additionally, the XGBoost analysis identified 13 important genes (Figure 5G). Finally, by integrating the results from all five ML approaches using a Venn diagram, we identified three overlapping genes, LINC-PINT, LRBA, and SYK, as the most robust and reproducible markers (Figure 5H).

Figure 5 Machine learning-based identification of key genes related to anal fistula-associated B cells. A: Box plot showing the expression levels of 19 candidate genes in B cells from anal fistula (AF) and control samples. Four genes were significantly upregulated in the AF group; B: Support vector machine model performance: Accuracy peaked (0.71) and error rate reached a minimum (0.29) when 13 variables were used; C: Coefficient profiles from least absolute shrinkage and selection operator (LASSO) regression analysis of the 19 genes across a range of penalty values (logλ); D: LASSO cross-validation plot identifying three key genes; E: Variable importance ranking of the top 15 genes derived from random forest analysis; F: Boruta algorithm identified 6 important genes; G: Extreme gradient boosting model ranking of variable importance across all input genes; H: Venn diagram showing the intersection of genes identified by five machine learning algorithms, with three common genes emerging as core candidates. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001. AF: Anal fistula; LASSO: Least absolute shrinkage and selection operator; RF: Random forest; XGBoost: Extreme gradient boosting; SVM: Support vector machine.

Immune infiltration analysis

To elucidate the link between gene expression and immune cell infiltration, ssGSEA algorithm was used to determine the abundance of immune cells in the AF samples. Subsequently, we examined associations between the expression of LINC-PINT, LRBA, and SYK and various immune cell types to identify distinct immunological patterns (Figure 6A). Additionally, we evaluated associations between these genes and immune-related functional pathways (Figure 6B). To investigate the biological roles of these genes more precisely, we implemented a single-gene GSEA by stratifying the disease groups according to the median expression level of each gene. In the LINC-PINT high-expression group, significant enrichment occurred in the TGF-β and MAPK signaling pathways (Figure 6C and D). In the LRBA high-expression group, the primary immunodeficiency and T-cell receptor signaling pathways were enriched (Figure 6E and F). For SYK, we detected cytokine-cytokine receptor interactions and NK cell-cell-mediated cytotoxicity pathways to be enriched (Figure 6G and H). These findings strongly indicated that LINC-PINT, LRBA, and SYK are critical for immune regulation and immune cell activation in AF and associated immune mechanisms.

Figure 6 Immune correlation and functional enrichment analysis of LINC-PINT, LRBA, and SYK. A: Correlation between the expression of LINC-PINT, LRBA, SYK and the abundance of various immune cells, based on single-sample gene set enrichment analysis (GSEA) scores in anal fistula samples; B: Correlation heatmap of LINC-PINT, LRBA, SYK with immune-related functions, including cytokine receptors, antigen processing, and co-stimulatory molecules; C and D: Single-gene GSEA enrichment analysis for LINC-PINT showing significant enrichment in TGF-β signaling pathway and MAPK signaling pathway; E and F: GSEA results for LRBA indicating enrichment in Primary immunodeficiency and T cell receptor signaling pathway; G and H: GSEA results for SYK showing significant enrichment in cytokine-cytokine receptor interaction and Natural killer cell mediated cytotoxicity. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001. KEGG: Kyoto Encyclopedia of Genes and Genomes.

Construction of the AF-associated clinical predictive model

We investigated the clinical diagnostic potential of hub genes LINC-PINT, LRBA, and SYK using transcriptomic data derived from AF samples. A box plot showed that all three genes exhibited differential expression between the AF and control groups, and LINC-PINT was significantly upregulated in the AF group (Figure 7A). To evaluate the predictive efficacy of each gene individually, we conducted a ROC analysis, which resulted in AUC values of 0.814 for LINC-PINT, 0.844 for SYK, and 0.778 for LRBA, all higher than the diagnostic threshold of 0.7 (Figure 7B-D). Subsequently, we incorporated these three genes into a clinical prediction model, which was represented as a nomogram (Figure 7E). This model showed strong predictive performance with an AUC of 0.858 after ROC validation (Figure 7F). Additionally, decision curve analysis demonstrated a strong concordance between the predicted and actual probabilities indicating that the model possessed excellent predictive accuracy for AF (Figure 7G).

Figure 7 Diagnostic performance and predictive model construction based on LINC-PINT, LRBA, and SYK. A: Box plot showing the expression levels of LINC-PINT, LRBA, and SYK in anal fistula (AF) vs control group. LINC-PINT showed significantly higher expression in the AF group; B-D: Receiver operating characteristic (ROC) curves evaluating the diagnostic performance of individual genes: LINC-PINT [area under the ROC curves (AUC)= 0.814], SYK (AUC = 0.844), and LRBA (AUC = 0.778); E: Nomogram constructed using LINC-PINT, LRBA, and SYK to predict AF risk; F: ROC curve evaluating the combined predictive model, showing strong discrimination ability; G: Decision curve analysis indicating favorable clinical utility of the nomogram-based prediction model. ^aP < 0.01; ^bP < 0.001. AF: Anal fistula; ROC: Receiver operating characteristic; AUC: Area under the receiver operating characteristic curve; DCA: Decision curve analysis.

Cell trajectory and pseudotime analysis

UMAP-based dimensionality reduction effectively illustrated expression patterns of LINC-PINT, LRBA, and SYK across various cell types in the AF samples (Figure 8A-C). To elucidate the dynamic expression trajectories of these key genes in B cells, we performed pseudotime analysis on B cells extracted from Seurat subsets using a monocle package. The pseudotime trajectory plot revealed continuous developmental progression of B cells, categorized into nine distinct cell states with pseudotime values reflecting a transition from early to late differentiation stages (Figure 8D). LINC-PINT exhibited elevated expression during early pseudotime stage followed by a gradual decline, suggesting its potential role in early B-cell activation and signaling (Figure 8E). LRBA expression showed an overall downward trend indicating that its functions were limited to the early and intermediate stages, with reduced effects during the later stages of activation (Figure 8F). Conversely, SYK expression increased during the late pseudotime stages suggesting its role in terminal B-cell differentiation and immune suppression (Figure 8G).

Figure 8 Cellular localization and pseudotime dynamics of key genes in B cells. A-C: Uniform manifold approximation and projection plots showing the expression patterns of LINC-PINT, LRBA, and SYK across all cell types in anal fistula samples; D: Monocle-based pseudotime trajectory analysis of B cells, identifying a continuous developmental path divided into nine states, with color gradients representing inferred pseudotime progression; E-G: Expression trends of LINC-PINT, LRBA, and SYK along the pseudotime axis. LINC-PINT shows high expression in early pseudotime, declining thereafter, suggesting a role in early B cell activation. LRBA shows a gradual decrease in expression, indicating a role in early-to-mid B cell stages. SYK expression increases at later pseudotime stages, implying involvement in terminal differentiation or immunoregulatory functions of B cells. UMAP: Uniform manifold approximation and projection.

Cell communication analysis using cell chat

Using single-cell transcriptomic data, we categorized B cells into two distinct subgroups, LINC-PINT⁺ and LINC-PINT^- based on the expression of the LINC-PINT gene. Cell chat analysis revealed that LINC-PINT⁺ B cells demonstrated significantly enhanced intercellular communication characterized by an increased interaction frequency and stronger communication signals, particularly non-B cell populations intercellular interaction (Figure 9A and B). We further visualized outgoing and incoming signaling patterns across various cell types. Notably, LINC-PINT⁺ B cells exhibited increased incoming signaling, especially within classical immune-related pathways such as MIF, CXCL, GALECTIN, and PROS (Figure 9C). Ligand-receptor analysis indicated significantly elevated signaling activities in LINC-PINT⁺ B cells involving pathways such as MIF, TGF-β, VEGF, NAMPT, and IL-6 (Figure 9D). The MIF-CD74 + CXCR4/CD44 axis suggested a potential role in inflammatory amplification and immune cell chemotaxis. The TGFB1-TGFBR1/2 and VEGFB-VEGFR1 pathways are involved in vascular remodeling and tissue repair. Additionally, axes such as IL-16-CD4, GRN-SORT1, and NAMP-TINSR indicate an extensive role of B cells for regulating T-cell reactions, macrophage activation, and metabolic signaling. Collectively, LINC-PINT⁺ B cells exhibit augmented intercellular communication capabilities and pronounced activation of the immune and inflammatory pathways. This suggested these cells had a higher immunoregulatory activity and functional importance compared to other cell subsets in the pathophysiology of AF and the associated inflammatory conditions.

Figure 9 CellChat analysis of LINC-PINT⁺ B cells in anal fistula. A: CellChat network plots illustrating the number and strength of intercellular communications among different cell types. LINC-PINT⁺ B cells exhibited more extensive and stronger interactions than LINC-PINT-B cells; B: Scatter plot of outgoing vs incoming signaling strength across all cell types. LINC-PINT⁺ B cells show higher overall communication strength; C: Heatmaps of outgoing and incoming signaling patterns across cell types. LINC-PINT⁺ B cells are highly active in classical immune-related pathways; D: The ligand-receptor interactions between different cell types and LRBA^+/- B cells.

B cells were further stratified into LRBA⁺ and LRBA^- subpopulations based on expression of the LRBA gene. Cell chat analysis demonstrated that compared to LRBA^- B cells, LRBA⁺ B cells showed significantly greater numbers of interactions and stronger overall communication with other cell types (Figure 10A and B). These findings indicate that LRBA expression might enhance the capacity of B cells to engage in immune or structural signaling crosstalk. Subsequent examination of signaling patterns across cell types revealed that LRBA⁺ B cells exhibited increased incoming signaling activity particularly enrichment in canonical inflammatory, immune-regulatory, and cell-adhesion pathways, and in pathways governing B-cell survival and differentiation such as MIF, CXCL, PTN, GALECTIN, and BAFF (Figure 10C). Further analysis of ligand-receptor pairs (Figure 10D) identified significant enrichment of the MIF-CD74 + CXCR4/CD44, VEGFB-VEGFR1, GRN-SORT1, and IL-16-CD4 signaling axes in LRBA⁺ B cells. These results highlighted the role of LRBA in the modulation of immune cell activation, vascular function, chemotaxis, and metabolic communication.

Figure 10  CellChat analysis of LRBA⁺ B cells in anal fistula. A: CellChat interaction network showing the number of intercellular interactions among all cell types. LRBA⁺ B cells exhibited a greater number of interactions compared to LRBA^- B cells; B: Communication strength analysis revealing that LRBA⁺ B cells possess higher overall interaction strength with other cell types; C: Heatmaps of outgoing and incoming signaling patterns across all cell types. LRBA⁺ B cells demonstrated enhanced outgoing signaling activity, particularly in classic immune and inflammatory pathways; D: The ligand-receptor interactions between different cell types and LRBA^+/- B cells.

To explore the immunological implications of SYK expression in B cells, we divided B-cell population into SYK⁺ and SYK^- subgroups and performed cell-cell communication analysis. The results revealed that SYK⁺ B cells engaged in a significantly higher number of intercellular interactions and exhibited enhanced overall communication strength compared to their SYK counterparts (Figure 11A and B). This suggested that SYK expression augmented the intercellular signaling capacity of B cells. Further analysis of outgoing signaling patterns demonstrated that SYK⁺ B cells exhibited markedly higher signal output activity particularly within canonical immune- and inflammation-related pathways, such as the MIF, GRN, and IL-16 signaling axes (Figure 11C). Ligand-receptor interaction analysis revealed a significant enrichment of critical signaling pathways in SYK⁺ B cells including MIF-CD74 + CXCR4/CD44, VEGFB-VEGFR1, and GRN-SORT1 (Figure 11D). Notably, the MIF axis may be involved in mediating inflammatory amplification and tissue remodeling, underscoring the regulatory role of SYK in IME dynamics.

Figure 11  CellChat analysis of SYK⁺ B cells in anal fistula. A: CellChat interaction network showing the number of interactions among different cell types. SYK⁺ B cells exhibit a markedly higher number of interactions compared to SYK^- B cells; B: Communication strength analysis reveals that SYK⁺ B cells possess significantly stronger interactions with surrounding immune and stromal cell populations; C: Outgoing and incoming signaling heatmaps demonstrate that SYK⁺ B cells show enhanced outgoing activity in key inflammatory and immune pathways; D: The ligand-receptor interactions between different cell types and SYK^+/- B cells.

Upregulated B cells hub genes derived from AF tissues

To study the role of B cells in AF, B cells were isolated from fistula tissues and healthy controls. Using flow cytometry, we identified B cells using a CD3^-CD19⁺ gating strategy to exclude T cells and analyzed gene expression with RT-qPCR (Figure 12). The baseline clinical characteristics of patients included in the RT-qPCR analysis were comparable between the two groups, with no statistically significant differences observed (Supplementary Table 2). Prior to quantification, the RT-qPCR assays were validated as follows: Melting curve analysis confirmed the specificity of amplification, showing a single sharp dissociation peak for each target gene without evidence of primer dimers or non-specific products (Supplementary Figure 2). Amplification plots demonstrated typical exponential growth phases across all samples, indicating high amplification efficiency (Supplementary Figure 3). Results showed that LINC-PINT, LRBA, and SYK mRNA levels were significantly higher in B cells from the fistula tissues indicating B-cell activation in the inflammatory environment.

Figure 12  Flow cytometric identification and gene expression analysis of B cells. A: Gating strategy for B cell isolation in healthy controls; B: Gating strategy for B cell isolation in samples of anal fistula. Cells were identified as CD19⁺ CD3^- to exclude T cells and ensure high purity of the B cell population; C: Reverse transcription-quantitative PCR analysis of LINC-PINT, LRBA, and SYK expression levels in the sorted CD19⁺ CD3^- B cells. ^aP < 0.01. AF: Anal fistula; HC: Healthy control.

DISCUSSION

AF is a chronic perianal inflammatory disease with a complex etiology that involves intricate interactions between immune and non-immune cell types. This complexity poses significant challenges for traditional research methodologies to accurately capture the highly heterogeneous nature of IME. ScRNA-seq allowed precise quantification of aberrantly infiltrated immune cell subsets under pathological conditions and aided identification of key molecular pathways within specific cellular subpopulations related to immune regulation, tissue repair, or inflammatory injury. Nevertheless, comprehensive analyses of B cells particularly those that employed scRNA-seq to explore B-cell-associated biomarkers for profiling immune cell infiltration and guiding precision therapy in AF are scarce.

Our study addressed this gap by utilizing scRNA-seq integrated with hdWGCNA to identify and characterize B-cell-associated subpopulations. Subsequently, a functional enrichment analysis of the genes within the identified co-expression modules was performed. To augment the sensitivity to gene expression levels and enhance the resolution of cellular heterogeneity, we combined bulk RNA-seq with single-cell transcriptomic data. This integrative approach enhanced data robustness and biological interpretability, facilitated a comprehensive understanding of gene expression dynamics, and enabled cell-level validation of large-scale transcriptomic findings to support the development of targeted therapeutic strategies[40]. Utilizing five ML algorithms, we identified three critical genes, LINC-PINT, LRBA, and SYK, which could be possible biomarkers for disease diagnosis and IME profiling. The diagnostic utility of these genes in AF was assessed through ROC curve analysis and a multivariable logistic regression-based predictive model which demonstrated excellent predictive accuracy for AF.

LINC-PINT is a member of the family of long noncoding RNAs (lncRNAs), which are RNA transcripts longer than 200 nucleotides that lack protein-coding potential[41]. Despite their lack of protein-coding capacity, lncRNAs are integral to gene regulation and function primarily in the nucleus. They communicate with RNA, DNA, and proteins to facilitate chromatin remodeling, transcriptional regulation, and post-transcriptional modifications[42,43]. Specifically, LINC-PINT has been implicated in the modulation of cellular migration and epithelial-mesenchymal transition in various cancer types[44]. Its expression is modulated by inflammatory signals, viral infections, and other cellular stimuli[45,46]. Importantly, the elevated expression of LINC-PINT in immune cells indicated its potential role in regulating IME, particularly in the context of inflammatory diseases.

The GSEA in our study also showed that the high expression of LINC-PINT was closely related to TGF-β and MAPK signaling pathways enrichment. In models of pancreatic ductal adenocarcinoma, LINC-PINT overexpression induced pronounced upregulation of TGF-β1 expression indicating that it might act as an upstream regulator of TGF-β[47]. Further, TGF-β is known to initiate both Smad-dependent and Smad-independent signaling pathways including the MAPK and PI3K/AKT pathways. These pathways facilitate the epithelial-to-mesenchymal transition of epithelial cells by inducing loss of cell polarity and intercellular junctions promoting the acquisition of mesenchymal migratory and invasive characteristics, and enhancing extracellular matrix deposition and fibrosis[48,49]. Collectively, these processes contribute to fistula progression and persistence. Furthermore, lncRNAs can function as competing endogenous RNAs (ceRNAs) that regulate microRNAs. In rheumatoid arthritis models, LINC-PINT acts as a ceRNA to sponge miR-155-5p, thereby alleviating its suppression of SOCS1. This, in turn, reduced the expression of inflammatory mediators such as IL-1β and matrix metalloproteinases, and inhibited the activation of the ERK/MAPK signaling pathways[50]. Accordingly, an elevated expression of LINC-PINT in patients with AF may not only enhance local fibrosis via upregulation of the TGF-β signaling cascade but also modulate MAPK pathway activity through ceRNA mechanism influencing B-cell activation and inflammatory responses. Collectively, our findings suggested that LINC-PINT might regulate chronic inflammation and fibrosis in AF.

LRBA is a signaling scaffold protein that orchestrates the assembly and localization of multiple signaling complexes in immune cells, serves as a key regulator of immune homeostasis, and protects against autoimmunity and intestinal inflammation. In this study, we found that LRBA expression in B-cell subpopulations was significantly elevated in the AF group compared with healthy controls. Previous studies showed that LRBA modulates the expression and recycling of CTLA4 in regulatory T cells (Tregs), thereby maintaining intestinal immune tolerance. LRBA deficiency results in Treg dysfunction and reduced CTLA4 expression collectively leading to the precipitation of intestinal inflammation, immune dysregulation, impaired mucosal barrier repair, and antimicrobial defense[51]. In addition, LRBA regulates BCR signaling cascades. The absence of LRBA leads to aberrant activation or phosphorylation of downstream BCR signaling molecules, including Btk, Plcγ2, IκBα, and p50, resulting in chronically activated B-cell phenotype characterized by reduced proliferative capacity and decreased survival[52]. Therefore, we speculate that LRBA upregulation in B-cell subpopulations in AF might represent a compensatory regulatory mechanism that supports B-cell functionality and contributes to the maintenance of local immune homeostasis.

This study found significant upregulation of SYK expression in B cells derived from individuals with AF. Immune infiltration analysis showed that the SYK-high group exhibited enrichment in two critical immune pathways, the cytokine-cytokine receptor interaction and NK-cell-mediated cytotoxicity, suggesting a highly activated local IME in AF. Research on inflammatory bowel disease (IBD) provide insights on the role of SYK in this context. As a pivotal tyrosine kinase involved in BCR signaling, SYK plays a crucial role in orchestrating B-cell activation, progression, and the production of proinflammatory cytokines. In the inflamed intestinal tissues of IBD patients both SYK expression and activity are significantly elevated, facilitating the secretion of inflammatory mediators such as TNF-α, IL-1β, and IL-6 from macrophages and other immune cells. This amplification of the inflammatory cascade contributes to disruption of the epithelial barrier[53-55]. In animal models, pharmacological inhibition of SYK effectively reduces colonic inflammation and restores barrier integrity, underscoring its potential as a therapeutic target[56]. Furthermore, SYK is integral to the function of NK cells as it regulates receptor-mediated signaling and cytotoxic responses. Although increased NK cell activity aids in the elimination of infected or abnormal cells, excessive NK cell activation can lead to heightened local tissue damage and inflammation[57,58]. In the context of AF, simultaneous overexpression of SYK and enrichment of NK-cell-mediated cytotoxic pathways indicate that innate and adaptive immunity assist in chronic inflammation and tissue destruction. The elevated expression of SYK in B cells suggests a locally activated, proinflammatory immune state similar to the SYK-driven inflammatory mechanisms observed in IBD. By amplifying cytokine signaling networks and enhancing NK cell cytotoxicity, SYK facilitates immune-mediated tissue damage within the AF microenvironment, highlighting its potential as a target for immunomodulatory therapy. Future studies should investigate the therapeutic potential of SYK inhibition in the management of AF-associated inflammation.

This study found LRBA and SYK were significantly upregulated in B-cell subpopulations from the AF group compared with healthy controls, and the results from RT-qPCR and scRNA-seq ere consistent. However, bulk transcriptomic data revealed that the overall LRBA and SYK expression was lower in the AF group than in the control group. This discrepancy in results could be attributable to the combined effects of complex cellular composition of fistula tissue and the intrinsic adaptive responses of individual cells. The results of scRNA-seq which focused on specific immune cell populations, suggested that under chronic inflammatory stimulation, infiltrating or tissue-resident B cells may adaptively upregulate LRBA and SYK. LRBA is essential for maintaining B-cell homeostasis as a critical regulator of intracellular trafficking, LRBA is essential for maintaining B-cell homeostasis. Its upregulation may represent a compensatory regulatory mechanism for preserving B-cell functionality and contributing to local immune homeostasis.

However, AF tissues analyzed in this study exhibited a significant increase in the number of infiltrating proinflammatory cells including neutrophils and activated T and NK cells, as well as non-immune stromal cells (such as fibroblasts). These proinflammatory and stromal cell populations might demonstrate low or negligible expression of LRBA and SYK. Consequently, although the mRNA levels of LRBA and SYK were elevated within individual B cells, the overall tissue-level expression appeared to be reduced owing to the dilution effect imposed by the high proportion of low- or non-expressing cell subsets. Taken together, the discrepancy between single-cell and bulk transcriptomic data partly reflects the complex immunopathology of AF. B cells may undergo adaptive regulation by upregulating LRBA, whereas increased SYK expression facilitates B-cell activation. However, the overall downregulation of LRBA and SYK at the tissue level reflects a pathological remodeling of the local IME, which is characterized by proinflammatory infiltrates and fibrotic stromal expansion. Future investigations integrating spatial transcriptomics are warranted for precise mapping of LRBA and SYK-high B cells within fistula lesions, to correlate their spatial distribution with severity of inflammation, and to elucidate their functional roles as potential therapeutic targets.

Reciprocal interactions between the endothelial, stromal, and immune compartments shape the AF microenvironment and govern the coupled trajectories of inflammation, angiogenesis, and fibrotic remodeling, yielding a B-cell-regulated network. Cell-cell communication analysis identified LINC-PINT⁺ B cells as central nodes with high incoming signaling activity. These cells preferentially integrated canonical immune programs, including MIF, CXCL, and GALECTIN, thereby consolidating inflammatory recruitment/retention, chemotactic positioning, and glycan-dependent adhesion remodeling within a single cellular node. Ligand-receptor mapping further delineated B-centered couplings—MIF/TGF-β (LINC-PINT⁺ B cell-fibroblast), VEGF/NAMPT (LINC-PINT⁺ B cell-endothelial cell), and MIF/IL-6 (LINC-PINT⁺ B cell-monocytic cell)—that assembled a positive-feedback architecture linking inflammatory drive to fibroblast activation, angiogenesis to metabolic provisioning and myeloid amplification, respectively. This framework provided a mechanistic rationale for persistent inflammation despite ongoing tissue repair in AF. A deeper mechanistic dissection of ligand-receptor crosstalk within the AF microenvironment will yield new insights for diagnosis and therapy.

This study has several limitations that should be acknowledged. First, our analyses focused primarily on B-cell-specific alterations and mechanisms. Although our CellChat analysis revealed active signaling between B-cell subpopulations and T/NK cells through axes such as IL-16-CD4 and MIF-CD74 + CXCR4/CD44, suggesting functional crosstalk that may collectively sustain the chronic inflammatory milieu, a more comprehensive investigation of T-cell subsets, including Tregs, and NK cell cytotoxic responses in the AF microenvironment will be an important focus of our future studies. Furthermore, integrating spatial transcriptomics would enable precise mapping of the co-localization patterns among B cells, T/NK cells, and stromal populations within fistula lesions, thereby deepening our understanding of the spatially organized immune networks underlying AF progression. Second, all transcriptomic analyses were performed exclusively on tissue specimens without paired peripheral blood samples. Consequently, whether the upregulation of LINC-PINT, LRBA, and SYK in AF-associated B cells reflects local tissue-specific inflammation, systemic immune dysregulation, or a combination of both remains to be elucidated. Future investigations incorporating paired tissue and peripheral blood specimens with longitudinal follow-up are warranted to resolve this question and to evaluate the potential of these hub genes as circulating biomarkers for non-invasive AF diagnosis. Third, validation in larger, multicenter, independent cohorts is required to establish the generalizability of these findings. Additionally, future studies should integrate detailed clinical phenotyping, encompassing fistula classification, disease duration, recurrence history, and comorbidity status, to evaluate associations between biomarker expression and clinicopathological features for improved patient stratification. Finally, the present findings are confined to the transcriptomic level; protein-level validation of LRBA and SYK, together with functional studies such as gene silencing in AF tissue-derived B cells, is warranted to establish causality beyond correlative observations.

CONCLUSION

Collectively, our study systematically identified cellular subpopulations specific to AF. By using ML and multiomics analysis, we elucidated the pivotal role of B-cell-associated genes in characterizing the immune cell infiltration of AF.

ACKNOWLEDGEMENTS

We sincerely thank all the patients for their tissue donation.