All datasets were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). ScRNA-seq data were obtained from two datasets: GSE214695 (Garrido-Trigo et al., 2023), which includes samples from 6 healthy individuals and 6 patients with active UC, and GSE231993 (Du et al., 2023), containing 4 healthy controls and 4 active UC cases. Microarray datasets were derived from three datasets: GSE87466 (21 healthy controls and 87 active UC samples) (Li et al., 2018), GSE75214 (11 healthy controls and 74 active UC samples) (Vancamelbeke et al., 2017), and GSE36807 (7 healthy controls and 15 active UC samples) (Montero-Meléndez et al., 2013). In patients with active UC, tissue samples were taken from inflamed regions of the colon, while in the control group, samples were obtained from the normal colonic tissue of healthy individuals. Table 1 shows the accession numbers, platforms, and other details of the datasets.

RAW264.7 cells were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). They were maintained in DMEM supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin. Cell cultures were incubated at 37 °C in a humidified atmosphere containing 5% CO₂. ELISA kit for IL-1β (ZC-37974, Zhuocai, China); ELISA kit for IL-10 (ZC-37962, Zhuocai, China); SPI1 antibody (55100-1-AP, Proteintech, China); iNOS antibody (22226-1-AP, Proteintech, China); Arg1 antibody (16001-1-AP, Proteintech, China); Tubulin antibody (80762-1-RR, Proteintech, China); Anti-rabbit IgG (H + L) (14780, Cell Signaling Technology, United States).

We conducted single-cell transcriptome analysis on the GSE214695 and GSE231993 datasets, utilizing the Seurat R package (version 5.1.0) for data processing. ScRNA-seq data underwent quality control filtering to retain cells with 300–6,000 detected genes, <10% mitochondrial reads, >1,000 UMIs, and total counts below the 97th percentile (Supplementary Figure S1). Data processing followed standard Seurat workflow including normalization, variable gene identification, and data scaling (Supplementary Figure S2). After initial dimensionality reduction via PCA, the Harmony algorithm was applied to mitigate batch effects (Supplementary Figures S3, S4). An SNN graph was constructed based on the top 40 harmony dimensions, and Louvain clustering at 0.4 resolution yielded 18 distinct cell clusters (Supplementary Figure S5). Finally, UMAP was used for two-dimensional visualization of cellular heterogeneity. Afterwards, we performed differential expression analysis on our Harmony-integrated scRNA-seq data using the FindAllMarkers function in Seurat, applying Wilcoxon rank sum test with stringent filtering criteria (logfc.threshold = 0.25, min.pct = 0.25, only.pos = TRUE). This approach enabled us to identify cluster-specific marker genes that facilitated accurate cell type annotation and provided valuable insights for subsequent functional characterization of cell subpopulations.

We leveraged pySCENIC methodology to elucidate GRN at single-cell resolution (Supplementary Figure S6). The analytical pipeline commenced with GRNboost2 algorithm implementation to swiftly detect potential regulatory interactions by examining transcription factor and target gene co-expression patterns within single-cell transcriptomic datasets (Van de Sande et al., 2020). We then incorporated the human promoter motif repository and used RcisTarget to identify binding sites within promoter regions, filtering for sequence-validated interactions to establish robust regulons. Regulatory influence across cellular populations was quantified via area under the curve cell-level (AUCell)-derived activity scores, complemented by PCA dimensionality reduction to offer alternative visualization of regulatory landscapes. The resultant cell-type-specific transcriptional networks illuminated distinctive molecular governance across diverse cellular states and developmental trajectories, providing crucial insights into heterogeneity and disease mechanisms. Moreover, we deployed cell-specific index (CSI) to characterize regulatory modules, with transcription factor relationships assessed through Pearson correlation and visualized.

To identify cell type-specific regulons, we first computed a binary indicator matrix representing the presence of each cell type across all cells. Using the regulon activity score matrix (rasMat), we then calculated the Regulon Specificity Score (RSS) for each regulon across different cell types. Specifically, RSS was defined as 1− Jensen-Shanno Divergence (JSD). 1− JSD between the regulon’s activity distribution and the idealized distribution of each cell type. A higher RSS indicates stronger specificity of a regulon to a particular cell type. The top-ranked regulons for each cell type were visualized based on their RSS values.

To investigate transcriptional regulation in UC, macrophage subsets were extracted from both UC and healthy control groups using preprocessed Seurat objects. Due to the sparsity of scRNA-seq data, a pseudobulk strategy was applied by aggregating gene expression profiles at the sample level, resulting in a gene-by-sample count matrix. Differential gene expression analysis was conducted using the DESeq2 package following median ratio normalization. Genes with at least 10 normalized counts in three or more samples per group were retained. Differentially expressed genes between UC and normal groups were identified based on a negative binomial generalized linear model. Gene set enrichment analysis (GSEA) was performed using the clusterProfiler R package (version 4.12.6) (Wu et al., 2021). Enriched pathways were then mapped onto the previously constructed transcription factor regulons.

To assess SPI1 activity in macrophages, we extracted its regulon target genes and calculated their average expression per cell using Seurat’s AddModuleScore (). This average expression served as an activity score reflecting SPI1 regulatory strength, which we visualized on UMAP and compared across macrophage subsets. Specifically, we assessed the expression activity of each SPI1 regulon target gene within macrophages and aggregated these measurements to estimate SPI1’s downstream regulatory impact.

To investigate the association between gene expression and clinical traits, we employed WGCNA to construct gene co-expression networks (Supplementary Figure S7). Initially, the WGCNA R package (version 1.72.5) was used to calculate pairwise gene correlations using Pearson’s correlation coefficient, resulting in a gene co-expression matrix (Langfelder and Horvath, 2008). To ensure the inclusion of biologically informative genes, we implemented the PVAC filtering strategy, which integrates variance, median absolute deviation (MAD), coefficient of variation (CV), and expression levels. Genes were retained for downstream analysis only if they exhibited high expression, fell within the top 75% in terms of variability, and had a minimum expression value greater than 1. Next, the pickSoftThreshold function was applied to determine the optimal soft threshold (β) for constructing a signed network, adhering to the scale-free topology criterion. This allowed us to construct a signed, weighted gene co-expression network that emphasizes positive correlations and conforms to a scale-free architecture, an important feature for identifying biologically meaningful modules associated with immune activation in UC. Modules were identified through hierarchical clustering and dynamic tree cutting, followed by module merging. Following module detection, we assessed the association between each module and UC by calculating the Pearson correlation between module eigengenes and the disease phenotype.

To identify WGCNA modules associated with SPI1 regulation, SPI1 target genes were first intersected with the set of expressed genes. The number of SPI1 targets within each module was then quantified. Subsequently, a hypergeometric test was conducted to evaluate the significance of target gene enrichment in each module. Modules exhibiting significant enrichment were deemed closely linked to SPI1 and were selected for downstream analyses. The module membership (MM) of each gene was evaluated, which is defined as the correlation between the gene and its module eigengene. The gene significance (GS) of a gene was also calculated, which is defined as the correlation with the UC phenotype. Genes with MM > 0.5 and GS > 0.2 were visualized via scatter plots, and those meeting both criteria were identified as candidate hub genes. Candidate hub genes identified from the module in one dataset were used to reconstruct a module eigengene in the other dataset via PCA. The association between the reconstructed eigengene and disease status was evaluated using Welch’s t-test. A significant difference in eigengene values supported the robustness of the association between the core module genes and UC.

To evaluate the reproducibility of co-expression modules, we conducted a module preservation analysis using WGCNA. Modules detected in GSE87466 were assessed for preservation in GSE75214 with 200 permutations and a signed network approach. Preservation strength was measured by the Zsummary statistic, where values above 10 indicate strong preservation, between 2 and 10 moderate preservation, and below 2 weak or no preservation. This analysis confirmed that key modules were consistently preserved across independent datasets, demonstrating their robustness.

Two microarray datasets, GSE87466 and GSE75214, were analyzed to identify differential expressed genes. A series of quality control procedures were applied to both datasets, including normalization using the normalizeBetweenArrays function, gene annotation and integration, and the detection of potential outliers through PCA (Supplementary Figures S8–S11). Following preprocessing, differentially expressed genes were identified using the limma R package (version 3.60.4), with the filtering criteria set to |logFC| ≥ 1 and a p-value <0.05 (Ritchie et al., 2015). Volcano plots were generated using the ggplot2 R package (version 3.5.1) to visualize the differential expressed genes. To investigate their biological relevance, GO enrichment analysis was conducted using the clusterProfiler package, focusing on processes associated with macrophage activation.

To identify potential downstream target genes regulated by transcription factors, LASSO, Recursive feature elimination with a random forest classifier (RFE-RF), and Support Vector Machine-Recursive Feature Elimination (SVM-RFE) were systematically applied. The training dataset was created by merging GSE75214 and GSE87466, with batch effects corrected using the sva R package (version 3.52.0) (Supplementary Figure S24). GSE36807 was designated as the independent validation set.

Given the high dimensionality of gene expression data, LASSO regression was first employed to perform variable selection and dimensionality reduction. Gene expression profiles were standardized, and LASSO logistic regression was conducted with 10-fold cross-validation to identify the optimal regularization parameter (λ) using the λ.1se criterion (Friedman et al., 2010). Genes with non-zero coefficients at this penalty level were considered candidate biomarkers with discriminative potential between UC and control samples. To ensure the robustness of selected features, stability selection was applied using the stabsel function, with a selection probability cutoff of 0.7 and PFER control (Hofner et al., 2015). Genes consistently selected across subsamples were considered stable and reliable predictors.

RFE-RF was implemented using the caret R package (version 6.0.94) (Liaw and Wiener, 2002), incorporating repeated 10-fold cross-validation to identify the optimal feature subset. Feature selection was guided by classification accuracy. A random forest model with 500 trees (ntree = 500) was then trained on the selected features, and variable importance was assessed via the Mean Decrease Gini index. To further confirm the robustness and relevance of the selected features, the Boruta algorithm was applied, with tentative attributes resolved using the TentativeRoughFix () function (Kursa et al., 2010).

In parallel, the SVM-RFE algorithm was employed to identify the most informative gene subset for classification, utilizing the caret package with a linear support vector machine classifier (svmLinear) implemented through the e1071 package (version 1.7.14) (Guyon et al., 2002). Prior to feature selection, gene expression data were row-wise standardized (Z-score normalization). The recursive feature elimination was performed with repeated 10-fold cross-validation (5 repeats) to ensure model stability and reduce overfitting. Various feature subset sizes were evaluated, and the optimal feature set was chosen based on the balance between classification accuracy and error rate.

To assess the diagnostic potential of the target genes, receiver operating characteristic (ROC) curves were generated using the pROC R package (version 1.18.5), and the corresponding area under the curve (AUC) values were computed (Robin et al., 2011). The GSE36087 dataset served as an external validation cohort. The expression levels of the target genes were analyzed in both UC and normal samples from the validation set, and ROC analysis was conducted to further confirm their discriminative ability.

RAW264.7 macrophages were subjected to M1-like polarization following two passages. Cells were treated with LPS (1 μg/mL) for 24 h to induce M1 polarization. After stimulation, cells were harvested and processed for subsequent assays. Throughout all experiments, untreated cells consistently served as blank controls to establish baseline values.

A small interfering RNA (siRNA) targeting SPI1 was designed and synthesized by Xuzhou General Biotechnology Co., Ltd. To evaluate the impact of SPI1 knockdown, cells were assigned to three groups: (1) control (untreated), (2) si-NC (transfected with scrambled, non-targeting siRNA as a negative control), and (3) si-SPI1 (transfected with siRNA specifically targeting SPI1).

Prior to transfection, cells were plated at 30%–40% confluence in antibiotic-free growth medium and cultured for approximately 24 h. Transient transfection was carried out using Lipofectamine reagent (Tiangen Biochemical Technology Co., Ltd.) following the manufacturer’s protocol.

All WB experiments were performed with three independent biological replicates (n = 3). Cells were lysed in RIPA buffer, and the resulting lysates were clarified by centrifugation. Protein concentrations were measured using the BCA assay and adjusted to ensure consistency among samples. Equal amounts of protein were loaded onto SDS-PAGE gels for separation and subsequently transferred to PVDF membranes. After blocking with 5% skim milk, membranes were incubated overnight with primary antibodies. Following washes, secondary antibodies were applied and incubated for 1 h. Protein signals were detected using an ECL substrate kit and quantified using ImageJ software.

Quantitative determination of inflammatory cytokines IL-1β and IL-10 was performed using ELISA methodology. All measurements were performed in triplicate using independent biological samples. Serial dilutions of reference standards were prepared in accordance with the manufacturer’s protocol to establish calibration curves. Absorbance measurements were obtained spectrophotometrically at 450 nm wavelength for all experimental samples, with blank wells serving as the reference baseline for instrument calibration. This standardized analytical approach enabled precise quantification of the target cytokines in the experimental specimens.

Flow cytometry experiments were conducted with at least three independent biological replicates. RAW 264.7 cells in the logarithmic growth phase were seeded into 6-well plates at a density of 1 × 10⁶ cells per well. After enzymatic digestion, cells were collected by centrifugation, resuspended in a fluorescent dye working solution, and incubated on ice for 15 min. Prior to analysis, cells were washed and resuspended in PBS. Flow cytometry was subsequently performed for detection.

Statistical analysis was performed using GraphPad 8.0 software. A Student’s t-test was used for two-group comparisons. For comparisons among more than two groups, a one-way ANOVA followed by a Tukey’s post hoc test was applied. Results are expressed as mean ± standard deviation (SD). Differences were considered statistically significant at P < 0.05 and highly significant at P < 0.01.

Although the two scRNA-seq datasets were derived from partially overlapping regions of the gut, we independently analyzed each dataset and observed a consistent increase in macrophage proportions in UC samples (Supplementary Figure S12). Given this consistent trend across both datasets, we integrated them for downstream analyses. This integration not only provides the increased cell numbers and diversity needed for regulatory network inference using pySCENIC, but also enables more comprehensive characterization of cell-type composition in the UC microenvironment.

Based on marker gene expression, the integrated scRNA-seq dataset was categorized into 14 distinct cell subpopulations (Figure 1A). Annotation of cell populations was performed using specific markers from the CellMarker database (Supplementary Figure S13). Analysis of cellular composition revealed distinct patterns between groups (Figure 1B). The control group was dominated by plasma cells (31.3%), followed by CD4 T cells (17.2%), CD8 T cells (13.3%), and naive B cells (13.0%), with lower proportions of fibroblasts (9.2%), epithelial cells (8.3%), and macrophages (2.1%). In contrast, UC samples showed significant shifts in cellular distribution. Immune cells expanded, with plasma cells increasing to 40.2%, CD4 T cells to 19.1%, and macrophages to 3.8%. Simultaneously, structural cells declined markedly, with fibroblasts decreasing to 2.9% and epithelial cells to 1.6%. Both CD8 T cells and naive B cells slightly decreased to 11.1% each, while all other cell types collectively comprised less than 11%.

Based on pySCENIC, we identified significantly active regulons. The activity of each regulon within individual cells was evaluated using the AUCell scoring algorithm, followed by UMAP dimensionality reduction clustering based on activity scores across all regulons (Figures 1C,D). Our analysis revealed notable differences in macrophages between UC and normal tissues. Based on the CSI among 401 regulons (TF-target gene pairs) identified by pySCENIC, cells were clustered into 9 modules (Supplementary Figure S14). The average activity of regulons within each module is illustrated in Supplementary Figure S15. Following PCA-based dimensionality reduction, a significant difference in macrophage enrichment between the UC and normal groups was observed (Figures 1E,F), suggesting that macrophages represent key differential cell types distinguishing UC from healthy controls.

To investigate the molecular mechanisms underlying altered macrophage function in UC, we conducted GSEA to identify signaling pathways enriched in macrophages. The analysis revealed significant enrichment of macrophage activation-associated pathway in UC compared to healthy controls (Figure 2A). To uncover potential upstream regulatory mechanisms driving these changes, we performed transcriptional regulatory network analysis. This revealed that SPI1 and its downstream targets comprised a substantial portion of the identified functional gene sets (Figure 2B). Furthermore, SPI1 regulon activity differed markedly between macrophages from UC patients and those from the control group (Figures 2C–E; Supplementary Table S1). We extracted SPI1 regulon target genes and calculated their average expression per cell using Seurat’s AddModuleScore () function. This average expression score served as a quantitative proxy for SPI1 regulatory strength. Visualization on UMAP embeddings showed heterogenous distribution of SPI1 regulon activity, with macrophages displaying significantly higher scores compared to other cell types (Figure 2F). We also calculated the average expression level of each SPI1 regulon gene within macrophages. Comparison between UC and control groups revealed that SPI1 target genes were significantly upregulated in UC macrophages, further supporting transcriptional activation of SPI1 (Figure 2G). These findings suggest that SPI1 is a key transcription factor involved in regulating macrophage activation in UC.

We applied the limma package to analyze two gene expression datasets, GSE87466 and GSE75214 (Figures 3A,B; Supplementary Tables S2, S3). Based on the identified differentially expressed genes, GO enrichment analysis was performed to investigate functional associations (Figures 3C,D). Both datasets revealed a significant enrichment of biological processes associated with macrophage activation.

To further explore the molecular mechanisms underlying UC, we applied WGCNA to the microarray gene expression profiles of both datasets. In GSE87466, one outlier sample (GSM2332144) was excluded based on hierarchical clustering (Supplementary Figure S16). A soft-threshold power of β = 14 was selected to construct a scale-free network (Figure 3F; Supplementary Figure S17). After module merging and refinement, 24 co-expression modules were identified (Figure 3H). Based on the enrichment analysis of SPI1 target genes, we identified the MEpaleturquoise module as significantly enriched and thus selected it for further investigation (Supplementary Figure S18; Supplementary Table S7). This module also exhibited a strong positive correlation with the UC phenotype (cor = 0.67, P = 2e−15), suggesting a potential association with macrophage activation in UC (Figure 3H). Genes within this module that met the criteria of MM > 0.5 and GS > 0.2 were defined as core genes (Figure 3J; Supplementary Table S4). In GSE75214, all samples clustered appropriately without any outliers, as shown by hierarchical clustering (Supplementary Figure S19).Similarly, in GSE75214, a soft-threshold of β = 12 was chosen (Figure 3E; Supplementary Figure S20), resulting in 13 modules after dynamic tree cutting and merging (Figure 3G). The MEgreen module was chosen due to its significant enrichment of SPI1 targets and moderate association with UC traits (cor = 0.58, P = 5e−9) (Supplementary Figure S21; Supplementary Table S6). Core genes from this module were similarly identified based on MM > 0.5 and GS > 0.2 (Figure 3I; Supplementary Table S5). These gene sets from both datasets were regarded as candidate targets potentially involved in macrophage-mediated pathogenesis of UC.

As shown in Supplementary Figure S22, most modules exhibited Zsummary >2, indicating moderate preservation across the two independent datasets. These findings demonstrate that key co-expression modules are preserved between GSE87466 and GSE75214, indicating cross-dataset consistency in gene co-expression patterns.

To evaluate the robustness of eigengene-trait associations, we performed cross-dataset validation between GSE75214 and GSE87466 through PCA. The association with disease status was tested using Welch’s t-test. The green module eigengene from GSE75214 showed a strong association with UC in GSE87466 (t = −13.96, P < 2.2e-16) (Figure 3K). Similarly, the paleturquoise module from GSE87466 was significantly associated with UC when projected into GSE75214 (t = −13.03, P = 1.42e-14) (Figure 3L). These results confirm that eigengene–trait relationships are consistent across datasets, supporting the stability of SPI1-enriched gene modules in UC.

We investigated the target genes regulated by SPI1 in macrophages and conducted an intersection analysis with two sets of core genes identified through WGCNA. This analysis revealed 272 overlapping core genes (Supplementary Table S8). To further refine the selection, LASSO regression identified 11 candidate genes with non-zero coefficients at the optimal regularization parameter (λ.1se) (Figures 4A,B). Stability selection confirmed 5 of the 11 LASSO-selected genes as high-confidence features (Supplementary Table S9). RFE-RF was employed to identify the optimal subset of genes for classification, with the best performance achieved using 10 features. This subset yielded an accuracy of 97.29% and a Cohen’s kappa of 0.9031 under 10-fold cross-validation repeated five times (Figure 4C; Supplementary Table S10). A final Random Forest model constructed using these 10 genes demonstrated strong classification performance and provided interpretable feature importance (Figure 4D; Supplementary Table S11). Notably, all 10 RFE-selected genes were also confirmed as important by the Boruta algorithm, further supporting their relevance (Supplementary Figure S23; Supplementary Table S12). Concurrently, SVM-RFE identified a 30-gene signature that yielded high classification accuracy (up to 98.7%) with minimal error and robust performance (Figure 4E; Supplementary Tables S13, S14). Finally, we selected the intersection of genes screened by the three machine learning methods and determined that the target genes were IRAK3, IL1RN, CD55 and PEA15 (Figure 4F).

To better evaluate the diagnostic performance of the four identified targets for UC, ROC curve analysis was conducted on both the training and validation datasets. The AUC values exceeded 0.90 in the training set and were above 0.80 in the validation set, demonstrating strong discriminatory power (Figures 4I,J). In two independent datasets (GSE87466 and GSE75214), SPI1, along with IRAK3, IL1RN, CD55, and PEA15 showed significantly higher expression in the UC group compared to the normal controls (P < 0.05, Figures 4G,H), suggesting that these genes may be upregulated in response to disease-associated regulatory mechanisms. Meanwhile, SPI1 expression showed a correlation with IRAK3, IL1RN, CD55 and PEA15 in human samples (Supplementary Figure S27), supporting a potential regulatory role of SPI1 in promoting downstream target gene expression.

To elucidate the molecular mechanism by which SPI1 regulates macrophage polarization, a series of in vitro experiments were conducted. Initially, M1 polarization of RAW264.7 mouse macrophages was successfully induced using LPS. WB analysis confirmed the induction, evidenced by a marked increase in the M1 marker iNOS and a corresponding decrease in the M2 marker Arg1 (P < 0.05, Figures 5A,C,D). Flow cytometry further supported this observation, showing a significant upregulation of CD86 in the LPS-treated group, while CD206 levels remained unchanged between groups, reinforcing the effective induction of M1 polarization (P < 0.05, Figures 5E–G). Subsequent WB analysis revealed that SPI1 expression was significantly elevated in the LPS group compared to the control (P < 0.05, Figures 5A,B). To investigate the functional role of SPI1, we employed siRNA to selectively silence SPI1 expression. The knockdown efficiency was validated through WB (P < 0.05, Figures 5H,I). WB analysis also showed that after SPI1 knockout, the level of the M1 marker iNOS was even higher than that induced by LPS, while the M2 marker Arg1 was decreased (P < 0.05, Figures 5J-L). ELISA analysis revealed that SPI1 knockdown significantly increased IL-1β secretion while reducing IL-10 levels following LPS stimulation (P < 0.05, Figures 5Q,R). Flow cytometry also found that knocking out SPI1 promoted M1 macrophage polarization (Figures 5M-P). These findings suggest that SPI1 plays a complex role in modulating macrophage polarization.