Gene expression and transcription factor regulatory network analysis identified four osteoblast subtypes

To understand the molecular features of human osteoblasts, we integrated cell clustering information obtained from gene expression profiles (by Seurat) and regulon activity profiles (by SCENIC). Three cell clusters of osteoblasts were identified based on systematic differences in their gene expression profiles (Fig. 2A). According to their order in the developmental trajectory (Fig. 2B) and their dynamic changes of the expression of BM-MSC and osteoblast-related marker genes (BM-MSC-related markers: LEPR [3], VCAM1 [4], osteoblast-related markers: RUNX2 [5], BGLAP [6]; Fig. 2C), we labeled the early-stage cell cluster with high expression of BM-MSCs-related markers as preosteoblasts, and the late-stage cell cluster with a high expression of osteoblast-related markers as mature osteoblasts. “Intermediate osteoblast” refers to the cell cluster in the middle stage of the developmental trajectory. SCENIC cell clustering results further identified two subtypes of preosteoblast [preosteoblast-S1 and preosteoblast-S2; Fig. 2D]. The regulon active scores of these two cell subtypes differed substantially from each other (Fig. 2E). These results reflect the heterogeneity of two cell subtypes of preosteoblasts, and thus, we identified a total of four cell clusters using integrated cell clustering information based on both gene expression and regulon activity.

Fig. 2

Single-cell clustering and regulatory network inference. A Gene expression-based cell clustering results shown using the Seurat clustering layout. Marked out three cell clusters. B Cell developmental trajectory inference based on gene expression. Marked three cell clusters. The upper-right trajectory plot indicates the direction of pseudotime. C Expression levels (log-normalized) of indicated genes in the three osteoblast clusters with respect to their pseudotime coordinates. The x-axis indicates the pseudotime, while the y-axis represents the log-normalized gene expression levels. Black lines depict the LOESS regression fit of the normalized expression values. D Gene expression and TF regulatory network conjoint clustering results shown using the SCENIC clustering layout. Marked out four cell clusters. E Regulons active heatmap in four cell clusters

Target gene function of active regulons in preosteoblast-S1

According to the transcription factor regulatory network activity heatmap (Fig. 2E), the 17 regulons demonstrating the highest active scores in preosteoblast-S1 were HES1, FOXC1, CEBPB, MAFF, FOSL2, CREM, EGR1, JUNB, NFIL3, ATF3, FOS, FOSB, CEBPD, JUND, BCLAF1, IRF1 and RUNX1. GSEA results showed that gene members in these regulons were mainly enriched in immunity and cell proliferation/differentiation-related function terms such as osteoblast differentiation (Fig. 3A), positive regulation of cell population proliferation, regulation of cell differentiation, regulation of cell population proliferation, regulation of cell activation (Additional file 2: Fig. S2A), inflammatory response (Fig. 3B), response to cytokine, leukocyte migration, cytokine receptor binding, and regulation of lymphocyte activation (Additional file 2: Fig. S2B). We used Cytoscape for network visualization and to mark these categories of gene function (e.g., immunity vs. cell proliferation/differentiation; Fig. 3C). The average expression levels of 113 gene members related to immunity and cell proliferation/differentiation showed a gradual downward trend in the four cell clusters as they matured (Fig. 3D).

Fig. 3

Active regulons in preosteoblast-S1. A GSEA analysis results in target genes of active regulons in preosteoblast-S1, “osteoblast differentiation” term. B GSEA analysis results in target genes of active regulons in preosteoblast-S1, terms of inflammatory response. C Network visualization of active regulons in preosteoblast-S1. Large green dots represent TF, small red dots represent immunity-related target genes, small blue dots represent cell proliferation/differentiation-related target genes, small yellow dots represent both immunity and cell proliferation/differentiation-related target genes; other target genes are shown as gray diamonds. D Violin plots showing average expression levels of gene members related to immunity and cell proliferation/differentiation. Dot size represents pseudotime for each cell from early (small) to late (large). *Adjusted p value < 0.05 (Kruskal–Wallis test)

Dynamic changes of CREM and FOSL2 regulon activity

The CREM regulon showed the highest active score (Fig. 2E) in preosteoblast-S1. There were 16 target genes in this regulon. Among these genes, ISG20, CXCL3, LIF, ABL2, IL6, and MTUS1 were related to immunity and/or cell proliferation/differentiation (Fig. 4A). Both CREM regulon activity and CREM gene expression showed specific high levels in preosteoblast-S1 (Fig. 4B, C). The average expression for target genes and expression levels for ISG20, CXCL3, LIF, ABL2, IL6, and MTUS1 were reduced gradually among the four cell clusters (preosteoblast-S1, preosteoblast-S2, intermediate osteoblast, and mature osteoblast, Fig. 4D, E). Target genes of the CREM regulon showed high expression levels at the early stage in pseudotime (Fig. 4F). The binding motif for most target genes of CREM is Taipale_cyt_meth_ATF4_GGATGACGTCATCC_eDBD (Fig. 4G, Table 1).

Fig. 4

Dynamic changes of CREM regulon activity. A CREM regulatory network. Red dots represent immunity-related target genes, blue dots represent cell proliferation/differentiation-related target genes, yellow dots represent both immunity and cell proliferation/differentiation-related target genes; other target genes are shown as gray dots. B CREM regulon activity, embedded on SCENIC clustering layout and colored by CREM regulon active levels. Gray dots represent completely inactive cells with zero active scores. Other dot colors represent active score in each cell from black (low) to blue (high). C Violin plots showing expression levels of CREM. *Adjusted p value < 0.05 (Kruskal–Wallis test). D Violin plots showing average expression levels of all target genes in CREM regulon. *Adjusted p value < 0.05 (Kruskal–Wallis test). E Immunity and cell proliferation/differentiation-related target gene expression in CREM regulon among four cell clusters. Dot color indicates the relative expression levels and the dot size shows the proportion of cells expressing the indicated genes. F Continuum of dynamic target gene expression in pseudotime of osteoblasts. Pixel color indicates the expression levels. Early-stage cells are list in on the left. G Motif taipale_cyt_meth_ATF4_GGATGACGTCATCC_eDBD of CREM. The relative sizes of the letters indicate the frequency of four bases in the sequences

Table 1 Motif information in CREM, FOSL2, FOXC2, RUNX2, and CREB3L1 regulons

FOSL2, which was found to be an upstream regulation gene of CREM, was the second most active regulon (Fig. 2E) in preosteoblast-S1. 14 of 39 FOSL2 target genes were related to immunity and/or cell proliferation/differentiation (Fig. 5A). Both FOSL2 regulon activity and FOSL2 gene expression showed specific high levels in preosteoblast-S1 (Fig. 5B, C). Average expression for FOSL2’s target genes decreased gradually over pseudotime (Fig. 5D). 14 immunity and/or cell proliferation/differentiation-related genes were highly expressed in preosteoblast-S1 (Fig. 5E). Target genes of the FOSL2 regulon also showed high expression levels at the early stage in pseudotime (Fig. 5F). The binding motif for most target genes of FOSL2 is cisbp_M3083 (Fig. 5G, Table 1). These collective results demonstrate that the CREM regulon, FOSL2 regulon, and their target genes related to immunity and/or cell proliferation/differentiation were highly active in the preosteoblast-S1 cluster.

Fig. 5

Dynamic changes of FOSL2 regulon activity. A FOSL2 regulatory network. Red dots represent immunity-related target genes, blue dots represent cell proliferation/differentiation-related target genes, yellow dots represent both immunity and cell proliferation/differentiation-related target genes; other target genes are shown as gray dots. B FOSL2 regulon activity, embedded on SCENIC clustering layout and colored by FOSL2 regulon active levels. Gray dots represent completely inactive cells with zero active scores. Other dot colors represent active score in each cell from black (low) to blue (high). C Violin plots showing expression levels of FOSL2. *Adjusted p value < 0.05 (Kruskal–Wallis test). D Violin plots showing average expression levels of all target genes in FOSL2 regulon. *Adjusted p value < 0.05 (Kruskal–Wallis test). E Immunity and cell proliferation/differentiation-related target genes expression in FOSL2 regulon among four cell clusters. Dot color indicates the relative expression levels and dot size shows the proportion of cells expressing the indicated genes. F Continuum of dynamic target gene expression in pseudotime of osteoblasts. Pixel color indicates expression levels. Early-stage cells are listed on the left. G Motif cisbp_M3083 of FOSL2. The relative sizes of the letters indicate the frequency of four bases in the sequences

To identify core target genes in the 17 highly active regulons we identified in preosteoblast-S1, we used immunity and cell proliferation/differentiation-related target genes to construct PPI networks (Fig. 6A, B). Interactions in these networks included correlation/regulation relationships or protein binding validated in Co-IP or other experiments [32]. Combined scores were used to value the confidence of interactions based on these evidences. High combined score interactions showed in wide edges were more valid than other interactions. Larger nodes with the higher network degrees represent the widely association of the corresponding proteins. Node color showed the gene expression ratio of preosteoblast S1 in comparison with mature osteoblast. Genes of the red/blue node were more highly expressed in preosteoblast S1 compared with other genes. MCODE extracted one core gene module in the immunity network with a MCODE score of 9.4 (Fig. 6C). Two core gene modules were extracted in the cell proliferation/differentiation network with MCODE scores of 5.4 and 4.7, respectively (Fig. 6D). Genes included in these core modules may widely influence the function of the whole network with more connection degrees. IL6 and TNFAIP3 were hub genes with high degrees in immunity network core module. IL6, CCL2, and PTGS2 were hub genes with high degrees in cell proliferation/differentiation network core modules. Only IL6 demonstrated the highest degree values in both immunity and cell proliferation/differentiation networks. In addition, IL6 also has many high combined score interactions and a high expression ratio of preosteoblast S1 compared with mature osteoblasts. These results comprehensively demonstrate that four hub genes, particularly IL6, had extensive interactions among core modules and might be critical for coordinating the function of target genes.

Fig. 6

PPI analysis of immunity and cell proliferation/differentiation-related target genes. A PPI network of immunity-related target genes. Dot size represents connection degrees from small (low) to large (high). Red nodes have the higher gene expression ratio of pre-osteoblast S1 in comparison with mature osteoblast. Edge width represents the combinded scores. B PPI network of cell proliferation/differentiation-related target genes. Dot size represents connection degrees from small (low) to large (high). Blue nodes have the higher gene expression ratio of pre-osteoblast S1 in comparison with mature osteoblast. Edge width represents the combinded scores. C Subnetwork in immunity-related PPI network screened by MCODE, MCODE score = 0.94. D Subnetwork in cell proliferation/differentiation-related PPI network screened by MCODE, MCODE score = 0.54 (left) and 0.47 (right). Dot color represents connection degrees in each subnetwork from blue (low) to red (high). Edge width represents the combined scores

Dynamic changes of MXD4 and KLF2 regulon activity

The MXD4 and KLF2 regulons showed higher active scores in preosteoblast-S2 (Fig. 2E). There were 29 target genes in the MXD4 regulon and 17 target genes in the KLF2 regulon (Fig. 7A, B). Although both MXD4 and KLF2 regulon activities showed specific high levels in preosteoblast-S2 (Fig. 7C), MXD4 gene expression was more specifically higher in the preosteoblast-S2 cell cluster (Fig. 7D). The average expression for target genes in the MXD4 regulon also showed the highest expression level in preosteoblast-S2 (Fig. 7E). Target gene clusters in the MXD4 and KLF2 regulons that showed high expression levels at the early stage in pseudotime are shown in Fig. 7F. The binding motif for most target genes of MXD4 and KLF2 is hocomoco__USF2_HUMAN.H11MO.0.A and transfac_pro__M07913 (Fig. 7G, H).

Fig. 7

Dynamic changes of MXD4 and KLF2 regulons activity. A MXD4 regulatory network. Target genes are shown as gray dots. B KLF2 regulatory network. Target genes are shown as gray dots. C MXD4 and KLF2 regulon activity, embedded on SCENIC clustering layout and colored by MXD4 or KLF2 regulon active levels. Gray dots represent completely inactive cells with zero active scores. Other dot colors represent active scores in each cell from black (low) to blue (high). D Violin plots showing expression levels of MXD4 (left) and KLF2 (right). *Adjusted p value < 0.05 (Kruskal–Wallis test). E Violin plots showing average expression levels of all target genes in MXD4 regulons. *Adjusted p value < 0.05 (Kruskal–Wallis test). F Continuum of dynamic MXD4 and KLF2 regulon’s target gene expression in pseudotime of osteoblasts. Pixel color indicates expression levels. Early-stage cells are list at left. G Motif hocomoco__USF2_HUMAN.H11MO.0.A of MXD4. H Motif transfac_pro__M07913 of KLF2. The relative sizes of the letters indicate the frequency of four bases in the sequences

Dynamic changes of FOXC2 and TAF7 regulon activity

Active scores for FOXC2 and TAF7 regulons were higher in intermediate osteoblasts than they were in the three other cell clusters (Fig. 2E). There were 26 target genes in the FOXC2 regulon and 11 target genes in the TAF7 regulon. EFNA1, SMAD7, SNAI1, HEY2, JAG1, SEMA7A and WNT4 in the FOXC2 regulon are genes related to MSC differentiation in the GO database (Fig. 8A, B). Both FOXC2 regulon activity and FOXC2 gene expression showed more specifically higher levels in intermediate osteoblasts than the other stages (Fig. 8C, D). Intermediate osteoblasts showed high expression levels of EFNA1, SMAD7, and SEMA7A (Fig. 8E) and average expression levels of all FOXC2 target genes (left plot in 8F). Target gene clusters in the FOXC2 and TAF7 regulons that showed high expression levels at the intermediate stage in pseudotime are shown in Fig. 8G. The binding motif for most target genes of FOXC2 and TAF7 is hocomoco__FOXL2_MOUSE.H11MO.0.C and dbcorrdb__TAF7__ENCSR000BLU_1__m1 (Fig. 8H, I, Table 1).

Fig. 8

Dynamic changes of FOXC2 and TAF7 regulon activity. A FOXC2 regulatory network. Blue dots represent MSC differentiation-related target genes; other target genes are shown as gray dots. B TAF7 regulatory network. Target genes are shown as gray dots. C FOXC2 and TAF7 regulon activity, embedded on SCENIC clustering layout and colored by FOXC2 or TAF7 regulon active levels. Gray dots represent completely inactive cells with zero active scores. Other dot colors represent active scores in each cell from black (low) to blue (high). D Violin plots showing expression levels of FOXC2 (left) and TAF7 (right). *Adjusted p value < 0.05 (Kruskal–Wallis test). E MSC differentiation-related target gene expression in FOXC2 regulon among four cell clusters. Dot color indicates relative expression levels and the dot size shows the proportion of cells expressing the indicated genes. F Violin plots showing average expression levels of all target genes in FOXC2 regulons. *Adjusted p value < 0.05 (Kruskal–Wallis test). G Continuum of dynamic FOXC2 and TAF7 regulon target gene expression in pseudotime of osteoblasts. Pixel color indicates expression levels. Early-stage cells are list at left. H Motif hocomoco__FOXL2_MOUSE.H11MO.0.C of FOXC2. I Motif dbcorrdb__TAF7__ENCSR000BLU_1__m1 of TAF7. The relative sizes of the letters indicate the frequency of four bases in the sequences

Activity of RUNX2 and CREB3L1 regulons increased as osteoblasts matured.

Active scores for RUNX2 and CREB3L1 regulons were higher in mature osteoblasts than they were in the three other cell clusters (Fig. 2E). Consequently, we focused on these regulons and performed GSEA to explore the function of their target genes. Although not significantly enriched, MEPE, PDGFC, MEF2C, and VDR in the RUNX2 regulon and PHOSPHO1, SPNS2, COL1A1 in the CREB3L1 regulon are genes related to skeletal system development in the GO database (Fig. 9A, B). RUNX2/CREB3L1 regulons were also relatively active in intermediate and mature osteoblasts (Fig. 9C). Expression levels for genes in the RUNX2 regulon rose gradually during osteoblast maturation, while CREB3L1 genes were significantly upregulated in mature osteoblasts (Fig. 9D). Intermediate and mature osteoblasts showed high expression levels of MEPE, PDGFC, MEF2C, and VDR (Fig. 9E) and average expression levels of all RUNX2 target genes (left plot in 9F). Moreover, high expression levels of PHOSPHO1, SPNS2, and COL1A1 were observed predominantly in mature osteoblasts (Fig. 9E). Average expression levels for all CREB3L1 target genes rose gradually during osteoblast differentiation, but expression levels increased dramatically in mature osteoblasts (right plot in 9F). Expression of most target genes in the RUNX2 regulon, including MEPE, PDGFC, MEF2C, and VDR, increased over pseudotime through the four cell subtype clusters (Fig. 9G). Despite specific target genes of CREB3L1 displaying a downward trend, genes related to skeletal system development (e.g., PHOSPHO1, SPNS2 and COL1A1) were gradually upregulated during osteoblast maturation (Fig. 9H). The binding motifs for most target genes of RUNX2 and CREB3L1 are swissregulon_hs_RUNX1.0.3.p2 and hocomoco_CR3L1_HUMAN.H11MO.0.D, respectively (Fig. 9I, J, Table 1). These results revealed that, as osteoblasts matured, there was an upward trend in RUNX2/CREB3L1 regulons and their target genes related to skeletal development.

Fig. 9

Dynamic changes of RUNX2 and CREB3L1 regulon activity. A RUNX2 regulatory network. Blue dots represent skeletal system development-related target genes, other target genes are shown as gray dots. B CREB3L1 regulatory network. Blue dots represent skeletal system development-related target genes, other target genes are shown as gray dots. C RUNX2 and CREB3L1 regulon activity, embedded on SCENIC clustering layout and colored by RUNX2 or CREB3L1 regulon active levels. Gray dots represent completely inactive cells with zero active scores. Other dot colors represent active score in each cell from black (low) to blue (high). D Violin plots showing expression levels of RUNX2 (left) and CREB3L1 (right). *Adjusted p value < 0.05 (Kruskal–Wallis test). E Skeletal system development-related target gene expression in RUNX2 and CREB3L1 regulon among four cell clusters. Dot color indicates relative expression level and the dot size shows the proportion of cells expressing the indicated genes. F Violin plots showing average expression levels of all target genes in RUNX2 (left) and CREB3L1 (right) regulons. *Adjusted p value < 0.05 (Kruskal–Wallis test). G Continuum of dynamic RUNX2 regulon target gene expression in pseudotime of osteoblasts. Pixel color indicates expression level. Early-stage cells are listed at left. H Continuum of dynamic RUNX2 regulon’s target gene expression in pseudotime of osteoblasts. Pixel color indicates the expression levels. Early-stage cells are listed at left. I Motif swissregulon_hs_RUNX1..3.p2 of RUNX2. J Motif hocomoco_CR3L1_HUMAN.H11MO.0.D of CREB3L1. The relative sizes of the letters indicate the frequency of four bases in the sequences

CSN analysis reflected connection changes between target gene pairs

Next, we analyzed important target gene connections in the four cell clusters at single-cell resolution. Gene connections in the CSN analysis are the inter-relationships (edges) among gene x and gene y in cell k that are assessed by statistic \(\hatx\left( k \right)\) (Eq. 1). We used target genes that are related to immunity, cell proliferation/differentiation, mesenchymal stem cell differentiation, and skeletal system development in the CREM, FOSL2, FOXC2, RUNX2, and CREB3L1 regulons to construct CSN for each cluster (Fig. 10A–C). Strong connections between 113 target genes in the CREM and FOSL2 regulons were established in preosteoblast-S1 (Fig. 10A). Preosteoblasts had the highest connectivity among 7 target genes in the FOXC2 regulon (Fig. 10B), while mature osteoblasts had the highest connectivity among 7 target genes in the RUNX2 and CREB3L1 regulons (Fig. 10C). Thus, immunity, cell proliferation/differentiation, and skeletal system development-related target genes in the CREM, FOSL2, RUNX2 and CREB3L1 regulons were only active in certain cell clusters. Their association network at single-cell resolution also has strong connections in the same cell types; however, strong connections of FOXC2’s target genes appeared before (in the early stage) the FOXC2 regulon attained the highest activation score (during the intermediate stage), which means that target genes in the FOXC2 regulon might be widely associated with each other in the early cell stage. These results suggest that FOXC2 may also play a role in the early stage of osteogenic differentiation mediated by target gene interactions.

Fig. 10

CSN construction among target genes. A CSN connections between immunity/cell proliferation/cell differentiation-related target genes in four cell clusters. Red arcs represent immunity-related target genes. Blue arcs represent cell proliferation/cell differentiation-related target genes. Yellow arcs represent immunity and cell proliferation/cell differentiation-related target genes. B CSN connections between skeletal system development-related target genes in four cell clusters. Purple arcs represent MSC differentiation-related target genes. C CSN connections between skeletal system development-related target genes in four cell clusters. Purple arcs represent skeletal system development-related target genes

Cell trajectory reconstruction reveals a new potential preosteoblast lineage

To further explore cell lineage from regulon activity aspect, we used regulons activity score matrix to reconstruct the regulon activity score-based cell developmental trajectory (Fig. 11A, B). We found that two preosteoblast subtypes, preosteoblast-S1 and preosteoblast-S2, were highly enriched in early cell stage and mature osteoblast was in the terminal stage (Fig. 11A, B). When we compared these results with the gene expression-based cell trajectory (Fig. 11C), we found that the uptrends of pseudotime values from preosteoblast, intermediate to mature osteoblast were coincident in both regulon score-based (Fig. 11A, B) and gene expression-based (Fig. 11C) cell trajectory results. Additionally, preosteoblast-S1 tended to form a distinct branch in the preosteoblast lineage (Fig. 11A, D), which was quite different in this trajectory structure. Branch heatmap of the branch point 2 in Fig. 11A, D (Fig. 11E) also showed a tendency toward upregulation of RUNX2 and CREB3L1 regulons in cell fate 1 (left branch after branch point 2 in Fig. 11D) and upregulation of CREM and FOSL2 regulons in cell fate 2 (right branch after branch point 2 in Fig. 11D). Although several of the 17 active regulons in preosteoblast-S1 showed a different activation tendency, this may be attributed to the mixed-cell composition in such cell fates (Fig. 11E). These results further strengthened the conclusion that preosteoblast-S1 was in a different cell state from the distinct branch of preosteoblast lineage compared with preosteoblast-S2, especially with regard to regulon activity.

Fig. 11

Cell trajectory reconstruction based on regulon activity. A Cell developmental trajectory inference based on regulon activity, four branch points in total. B The direction of pseudotime in trajectory plot for Fig. 9A. C Cell developmental trajectory branch plot based on gene expression. The upper-right trajectory plot indicates the direction of pseudotime. D Cell lineage relationships in Fig. 9A. E Continuum of dynamic target gene expression around branch point 2 in Fig. 9A, D. Cell fate 1 was correlated to the up branch after branch point 2 in Fig. 9A (left branch after branch point 2 in Fig. 9D). Cell fate 2 was correlated to the down branch after branch point 2 in Fig. 9A (right branch after branch point 2 in Fig. 9D). Pixel color indicates the expression levels. All target genes are clustered in 2 groups based on their expression pattern. F NDM of CREM, FOSL2, FOXC2, RUNX2, and CREB3L1 TFs. Stars indicate the significance levels of the NDM difference from any other cell clusters (adjusted p value < 0.05, Wilcoxon rank-sum test)

To explore gene connections at the whole transcriptome level, we calculated network degrees of five highlighted TFs. CREM and FOSL2 demonstrated the highest NDM values in preosteoblast-S1, FOXC2 demonstrated the highest NDM values in intermediate osteoblasts, while the NDM values of RUNX2 and CREB3L1 approached their peaks in mature osteoblasts (Fig. 11F). These results indicate that more interactions exist between these TFs and other genes in the corresponding cell subtypes. These are also the same trends identified for the activity of the corresponding regulons.

We also used diffusion mapping [28] to further confirm our osteogenic differentiation trajectory results based on gene expression and TF regulation aspects. Pseudotime order of preosteoblasts (preosteoblast-S1, preosteoblast-S2), intermediate osteoblasts, and mature osteoblasts in diffusion map analysis results were consistent with our earlier Monocle-based analysis (Fig. 12A–D). Compared with preosteoblast-S2, diffusion map analysis showed that preosteoblast-S1 cells were more concentrated in the DC2 dimension (Fig. 12A) and the early stage of pseudotime (< 1500, Fig. 12B) in gene expression-based trajectory. Preosteoblast-S1 also tended to form a distinctive branch in preosteoblast lineage in TF regulation-based trajectory (Fig. 12C, D). Thus, preosteoblast-S1 performed differently than preosteoblast-S2, even with the independent approach using gene expression-based “traditional” pseudotime analysis. These differences were more apparent in TF regulation-based trajectories, which was consistent with the results from the Monocle analysis. These results further validated the independence of preosteoblast-S1. CREM, FOSL2, FOXC2, RUNX2, and CREB3L1 regulon coincident active tendencies were also confirmed in monocle (Fig. 12E) and diffusion map trajectory heatmap (Fig. 12F). Thus, CREM and FOSL2 regulons were highly active in the early cell stage (preosteoblast-S1), FOXC2 regulon was highly active in the intermediate cell stage, and RUNX2 and CREB3L1 regulons were highly active in the late cell stage (mature osteoblast).

Fig. 12

Validation of regulon activity tendency. A Diffusion map trajectory inference based on gene expression. The upper-right trajectory plot indicates the direction of pseudotime. B Cell distribution based on the pseudotime coordinates (Fig. 10A). The x-axis is the pseudotime and the y-axis represents the osteoblast subtypes. C Diffusion map trajectory inference based on regulon activity. The upper-right trajectory plot indicates the direction of pseudotime. D Cell distribution based on the pseudotime coordinates in (Fig. 10C). The x-axis is pseudotime and the y-axis represents osteoblast subtypes. E Monocle trajectory heatmap of CREM, FOSL2, FOXC2, RUNX2, and CREB3L1 regulons. F Diffusion map trajectory heatmap of CREM, FOSL2, FOXC2, RUNX2, and CREB3L1 regulons. G Immunity and cell proliferation/differentiation-related target genes’ average expression in CREM regulon. H Immunity and cell proliferation/differentiation-related target genes’ average expression in FOSL2 regulon. I MSC differentiation-related target genes’ average expression in FOXC2 regulon. J Average expression of all target genes