We obtained three UC single-cell RNA-sequencing datasets [GSE116222; GSE114374; GSE95459] and one UC RNA bulk dataset [GSE165512] from GEO (http://www.ncbi.nlm.nih.gov/geo). Following standardization, samples deficient in clinical information were discarded (Guo et al., 2023). The following are nine samples in GSE116222 (3 Healthy, three adjacent non-inflamed areas of UC, and three inflamed areas of UC samples), four samples in GSE114374 (2 Healthy, 2 UC samples), and 10 samples in GSE95459 (5 Healthy, 5 UC samples). Additionally, 86 samples were obtained from GSE165512 (46 Healthy and 40 UC samples). Dataset GSE95459 is specifically designed for analyzing epithelial cell heterogeneity and gene expression in UC. All analyses related to immune cell populations (e.g., T cells, B cells) and IRGs (e.g., CD28) were conducted using the GSE116222 and GSE114374 datasets, which contain full-thickness mucosal biopsies encompassing diverse cell types. Table 1 displays the features of the four datasets. The clinical information of GSE116222 samples is provided in Supplementary Table S1.

We retained cells that expressed above 200 but no more than 4,000 genes. Meanwhile, a cutoff value of 5% of mitochondrial genes was established for further filtering. To produce cell clusters that could be shown and annotated using the t-SNE diagram, the number of principal components (PCs) was increased to 20 after 1,500 hypervariable genes were identified for study. The top 10 unique expression genes in each cluster were then selected using the “FindAllMarkers” function of the Seurat R package. Then a total of 11 clusters were found (Supplementary Table S2).

To automatically label our single-cell RNA-seq data, we used the R package “SingleR”. Between the expression profiles of each cell and those of the reference sample, Spearman’s correlation was calculated. According to a prior study, cell type identification is based on differentially expressed genes (DEGs) in each cluster with manual verification (Zhang et al., 2019). We applied this strategy repeatedly for each label, annotating the cell with the label that received the highest score. We found the genes that were differentially expressed between UC and normal cells using the “FindMarkers” approach. A list of each major indication that separates UC from normal cells may be seen in Supplementary Table S3. The reference marker gene list used for manual annotation is provided in Supplementary Table S4.

After each cell had been annotated, the inflammatory cell objects were extracted, and cells with mean expression >0.1 and dispersion empirical >1 * dispersion fit for the next stage of the pseudotime study (Xiong et al., 2022). We then used the “DDRTree” technique to decrease the dimension of the cells, using the reduceDimension function to identify the type of cell differentiation state. Finally, we used the plot_cell_trajectory function to depict the cell differentiation trajectory.

We used the limma package (version 3.8) to identify DEGs on the raw data of GSE165512 (Ritchie et al., 2015; Tang et al., 2022). Significantly differentially expressed genes (DEGs) are those with an absolute |log₂FC| >1 and adjusted p-value <0.05. With the help of the ggplot2 program, volcano and heatmap plots were created (version 3.6.3). After that, we carried out GO and KEGG analyses as well as ranking analyses of the DEGs. Table 1 provides a comprehensive list of all datasets used in this investigation.

On GSE165512, we ran a Weighted Gene Co-Expression Network Analysis (WGCNA) analysis (Guo et al., 2024) and then used a subset of genes with an expression standard deviation greater than 0 for additional analysis, removing outlier data. The data were separated into distinct modules by selecting an optimum soft threshold (Supplementary Table S5) and concurrently identifying the modules that were most closely related to UC.

The R program ConsensusClusterPlus was used to group the inflammatory cells into various categories. When the clustered index “k” increases from two to 9, choose an appropriate value of k that maximizes the difference between clusters. Then, we draw a box diagram with the R package ggplot2 and the R package reshape2. To explore whether there is a correlation between genes in branches and clusters. Using the CIBERSORT algorithm, we summarized the distribution of immune cell subsets in inflammation samples in UC. In addition, we further performed an analysis of immune cells and IRGs among different types.

For GO and KEGG analysis, the DEGs for the T cluster in scRNA sequencing data and DEGs in bulk sequencing data were separately imported into xiantao, an online bioinformatic analysis platform (https://www.xiantao.love/) (Zhou et al., 2022). Based on p-value ranking, the top 10 paths were chosen.

Protein-protein interaction (PPI) network analysis was performed using STRING (https://string-db.org/).

Colitis was induced in male C57BL/6 mice by treatment with 3% dextran sulfate sodium (DSS) for 7 days. Successful modeling was confirmed when the disease activity index (DAI) score exceeded 3, indicating weight loss, diarrhea, and bloody stools. Colon tissues were collected from DSS-induced colitis model mice (n = 6) and healthy control mice (n = 6), and total RNA was extracted using Trizol reagent (Thermo Fisher Scientific, MA, USA). 5 μg of total RNA was reverse-transcribed into cDNA following the manufacturer’s instructions. CD28 mRNA expression was measured using SYBR Green Master Mix (Promega, WA, USA) and the CFX 96 Real-Time PCR System (Bio-Rad Laboratories, CA, USA). Primer sequences are listed in Table 2 and were synthesized by Wuhan Service Biotech. GAPDH served as an internal reference, and relative expression levels were calculated using the 2 (−ΔΔCt) method.

Colon tissue slices from the colitis model mice (n = 6) and the healthy control mice (n = 6) were collected. After sacrificing the mice, tissue was washed with ice-cold phosphate-buffered saline (PBS), fixed with 4% formaldehyde, and embedded in paraffin. Tissue sections (4 μm thick) were stained with anti-CD28 antibody (clone number, EPR22076, Abcam) using immunohistochemistry (Fujii et al., 2021). Two researchers independently evaluated each section. Sections from each mouse were randomly examined using four different visual fields to calculate the average optical density value for each tissue section.

RNA-Seq differential expression was analyzed using limma with Benjamini–Hochberg FDR correction. WGCNA used Pearson correlation for module-clinical trait associations. IRGs and DEG intersections were analyzed with GO and KEGG enrichment (P < 0.01). Single-cell analysis used Seurat and SingleR, and immune cell infiltration was assessed with CIBERSORT. Group differences were tested with an unpaired t-test (P < 0.05).

An exercise in quality control was done on the single-cell dataset. To ensure the quality of the cell samples utilized in the research, as indicated in Figure 1A, we excluded specific cells with <200 genes, >4,000 genes, and controlled the fraction of mitochondrial genes. With a correlation value of 0.71, the nCount RNA, which indicates the number of distinct molecular identifiers, is inversely linked with the percentage of mitochondrial genes (Figure 1B). With a correlation coefficient of 0.89, the number of genes represented by nFeature RNA and nCount RNA is positively associated (Figure 1C). Thereafter, we identified 1,500 genes with high variability, all of which are indicated in red, and marked the 10 most important genes (Figures 1D,E). TPSB2 and TPSAB1 are the top two HVGs, which are neutral proteases present in mast cells. They are elevated in the colon of inflammatory bowel disease patients and may be involved in the innate immune response (Cenac et al., 2002). TPSB2 and TPSAB1 fold changes and significance are shown in Supplementary Table S6.

As shown in Figure 2A, we first performed PCA dimensionality reduction analysis to obtain the genes associated with each PC. Then we calculated the coordinates of each cell in PC-1 and PC-2 by PC correlation coefficient, and labeled them (Figure 2B). Using JackStrawPlot, PCA detected all 20 PCs with a p-value <0.05 (Figure 2C), and then the 20 PCs were analyzed by TSNE dimensionality reduction. The genes highly expressed in each PC were marked in yellow in the heat map (Figure 2D).

We found the genes that were significantly different between UC and normal cells using the “FindMarkers” approach. A previous study suggested that these clusters could be connected to known cell lineages using marker genes (Sinha et al., 2018). Visualization of 11 clusters using the t-SNE analysis (Figure 3A). T cells and B cells cluster revealed an increased percentage in the UC group (Figure 3B), and the T cells cluster was the main focus of the analysis that followed. The dot plot (Figure 3C) and violin plot (Figure 3E) display the expression of T cell type marker genes. In addition, we mapped t-SNE plots of T cell type marker genes in all clusters (Figure 3D). The DEGs in each cell cluster were heatmapped, and the top 10 DEGs were marked with a yellow sticker (Figure 3F).

We used simulation to analyze the cell trajectory differentiation of UC cells. We found that the darker the blue, the earlier the cell differentiation, and the lightest blue denoting the most recently differentiated cells (Figure 4A), indicating that inflammatory cells gradually differentiate from left to right. All inflammatory cells had six distinct developed states, each designated with a different hue, as seen in Figure 4B, with the red type being the earliest differentiated type. After that, we looked into how various cell clusters differentiate (Figure 4C). All cells were analyzed, including B cells, T cells, monocytes, and epithelial cells (Figure 4D). We discovered that T cells and B cells differentiated later than normal epithelial cells.

The bulk RNA sequencing dataset GSE165512, which includes 40 UC patients and 46 healthy controls, was used to examine the gene expression characteristics in UC. After the null value is eliminated, there are 28,864 IDs remaining. DEGs with |log₂FC| >1 and adjusted p-values <0.05 were chosen (Supplementary Table S7). Following that, 947 upregulated and 4,633 downregulated DEGs were kept (Figure 5A). The difference multiples and adjusted p-values were obtained after the difference analysis. The difference ranking chart was made by the size of the final difference multiples data to show the results of difference analysis (Figure 5B). Relative consistency was seen within groups in the heatmap of the top 30 upregulated and top 30 downregulated DEGs (Figure 5C). Next, we used the GO database to obtain the association of differentially expressed genes at three levels: biological process (BP), cellular component (CC), and molecular function (MF). Interestingly, the top two terms of BP, CC, and MF of DEGs in UC were mainly focused on immune response (Figure 5D). The KEGG pathway was enriched for Bile secretion, Fat digestion and absorption, Protein digestion and absorption, Metabolism of xenobiotics by cytochrome P450, Chemical carcinogenesis, Neuroactive ligand-receptor interaction, Steroid hormone biosynthesis, Drug metabolism cytochrome P450, ABC transporters, Serotonergic synapse (Figures 5E,F).

First, we discovered 993 genes using WGCNA on GSE165512. Second, the cluster analysis was carried out using the “flashClust” utility package; the samples were grouped into four statistically different groups. The clusters with a small number were removed, and the other three clusters were used for further analysis (Figure 6A). The “choose Soft Threshold” function of the “WGCNA” package was then used to filter the power parameter range of 1–20. For the purpose of building a scale-free network, we chose a power of b = 15 as the soft threshold (Figure 6B). We set the threshold to 0.3 (Figure 6C) and the minimum number of modules to 30 to combine related modules in cluster 3. Ten modules encompassing genes with comparable co-expression characteristics were created (Figure 6D). Multiple modules were related to UC, with the green module being the most important and encompassing 69 genes, as shown by module-trait association analyses (Figure 6E). The green module showed a significant association with both healthy samples and UC (COR = 0.64, P < 0.001). In the green module, we analyzed the top 10 WGCNA-hub genes (Figure 6F).

To explore the heterogeneity of UC, we clustered inflammatory cells in the published dataset GSE165512 using a previously developed consensus unsupervised clustering technique. When K = 3, there is a flattest middle portion of the CDF curve in the consensus matrix (Figure 7A). Additionally, we discovered that when K = 3 was chosen for the consensus clustering analysis, the interference between subgroups could be minimized. Therefore, the analysis defined three clusters with the most robust classification (Figure 7B). Next, we performed expression analysis of the typing genes. We aimed to observe the correlation between the genes in the branches of the pseudotime analysis in the single-cell analysis and the clusters in the genotyping. We found that genes upregulated in branch 1 were highly expressed in cluster 1, and genes downregulated in branch 1 were also under-expressed in cluster 1. A similar pattern was found in branch 5 and cluster 3 (Figure 7C). We initially compiled the findings from UC samples using the CIBERSORT method. The correlation heatmap of the 22 immune cell subpopulations in inflammatory samples shows this (Figure 7D). By immune cell analysis, we found that there were significant differences in naive B cells (P < 0.001), Plasma cells (P < 0.05), M2 Macrophages (P < 0.05), and resting Mast cells (P < 0.01) between different types of clusters (Figure 7E). In addition, we further performed immune checkpoint analysis to observe whether there were differences in immune checkpoint-related genes among different types. We found differences in these immune checkpoints, such as TNFRSF9, PVR, PTPRC, PDCD1, ICOS, CTLA4, CD8A, CD80, CD40LG, CD40, CD28, among which PVR, ICOS, CD28 was the most significant difference (P < 0.001) (Figure 7F).

To study the transcriptionally regulated activity of IRGs, we downloaded a list of human transcription factors from Human TFDB (http://bioinfo.life.hust.edu.cn/Human TFDB/#!/) (Hu et al., 2019). We screened 84 TFs in the UC T cluster in single-cell data and 222 TFs in the UC bulk data. Furthermore, the UC T cluster and the UC intestinal tissue both expressed 11 shared TFs (Figure 8A). With the default settings (only.pos = FALSE, min. pct = 0.25, and logfc. threshold = 0.5), the FindAllMarkers function was run. Across all clusters, 7,723 different marker genes were discovered. There were 1,416 marker genes found for the T cluster in total, however only 84 of these genes corresponded to TFs. There were only 11 overlapping TFs between the T marker genes, DEGs in the UC intestinal tract, and the TF database. The expression of 11 TFs in the data set GSE165512 and the single-cell T cluster is shown in Figures 8B,C. 5 TFs were upregulated in the UC intestinal, and 3 TFs were raised in T cells. As a hub gene for the transcriptional control of IRGs, the PPI network predicts that HNF4A may play a significant role (Figure 8D). The names of the 11 co-expressed TFs and their expression trends in UC and healthy samples are provided in Supplementary Table S8.

By using qRT-PCR, we were able to determine the biomarkers’ expression levels. CD28 showed a significant upregulation in DSS-induced model mice colonic tissues (P < 0.05) (Figure 9A). After that, the expression level of CD28 was detected by immunohistochemistry. CD28 was highly expressed in colonic tissues of DSS-induced colitis mice (P < 0.001) (Figures 9B,C).