Skip to main content


We’d like to understand how you use our websites in order to improve them. Register your interest.

Time series gene expression profiling and temporal regulatory pathway analysis of BMP6 induced osteoblast differentiation and mineralization



BMP6 mediated osteoblast differentiation plays a key role in skeletal development and bone disease. Unfortunately, the signaling pathways regulated by BMP6 are largely uncharacterized due to both a lack of data and the complexity of the response.


To better characterize the signaling pathways responsive to BMP6, we conducted a time series microarray study to track BMP6 induced osteoblast differentiation and mineralization. These temporal data were analyzed using a customized gene set analysis approach to identify temporally coherent sets of genes that act downstream of BMP6. Our analysis identified BMP6 regulation of previously reported pathways, such as the TGF-beta pathway. We also identified previously unknown connections between BMP6 and pathways such as Notch signaling and the MYB and BAF57 regulatory modules. In addition, we identify a super-network of pathways that are sequentially activated following BMP6 induction.


In this work, we carried out a microarray-based temporal regulatory pathway analysis of BMP6 induced osteoblast differentiation and mineralization using GAGE method. This novel temporal analysis is more informative and powerful than the classical static pathway analysis in that: (1) it captures the interconnections between signaling pathways or functional modules and demonstrates the even higher level organization of molecular biological systems; (2) it describes the temporal perturbation patterns of each pathway or module and their dynamic roles in osteoblast differentiation. The same set of experimental and computational strategies employed in our work could be useful for studying other complex biological processes.


Osteoblasts are responsible for bone matrix production and mineralization [1]. In concert with osteoclasts, osteoblasts coordinate bone remodeling, a physiologic process by which bone mass is maintained [1]. Osteoblasts arise from mesenchymal stem cells (MSC) [1]. However, the mechanism controlling this differentiation process is not well understood.

MSC differentiation to osteoblasts can be induced both in vivo and in vitro by soluble factors in the bone morphogenetic protein family (BMP) [2, 3]. Among the BMPs, BMP2, 4, 6 and 7 are the best characterized osetogenic factors [4]. Our previous work [5] has shown that: (1) human MSC produce BMP6 in defined, serum-free conditions, (2) BMP6 is up-regulated under mild osteogenic stimulus (dexamethasone), (3) exogenous BMP6 potently induces osteoblast differentiation, but responses to BMP2, 4, or 7 are inconsistent and require higher doses, (4) exogenous BMP6 induces the expression or up-regulation of a set of osteoblast-related genes in human MSC, and (5) 24 hour treatment with BMP6 induce high levels of osteoblast gene expression and cause mineralization. These results established BMP6 as an endogenous regulator of human osteoblast differentiation [5].

Unfortunately, BMP6 signaling largely remains uncharacterized. Significant efforts, particularly a series of high throughput microarray studies [612] have been undertaken to uncover BMP (including BMP6) responsive genes, transcriptional programs, and their roles in osteoblast development. However, an integrated understanding of the regulatory mechanisms for osteoblast differentiation and mineralization has not been achieved. In particular, two remaining questions include: (1) what pathways and gene groups are responsible for MSC differentiation to osteoblasts in response to BMP6 stimulation? (2) what is the sequence of pathway activation events in response to BMP6 during the osteogenic induction?

To answer these two questions, we conducted a time series microarray study on BMP6 osteogenic induction and a comprehensive pathway analysis on the temporal data. We met special challenges at both experiment and data analysis levels.

At experiment level, we considered two major issues. First, what rate to sample the time series? The time intervals should be short enough to capture the dynamics and continuity, but long enough to show phenotypically significant changes. We determined time intervals based on the minimal BMP6 treatment durations needed for significant phenotypic changes, including the expression of osteoblast markers and formation of mineralized extracellular matrix (Figure 1). Second, given that osteoblast induction is a temporal process with accumulative effects and gradual changes, how to dissect out the net effect of BMP6 at different phenotypic stages? We employed a BMP6 addition-and-withdrawal scheme (Figure 1).

Figure 1

Design for the microarray study on BMP6 induced osteoblast differentiation. Human MSC cells were pre-cultured for 4 days and subsequently treated with BMP6 for 0 hours, 8 hours, 24 hours, and 96 hours. These four time points correspond to four phenotypic groups, of control, preosteoblast (no mineralization), (sub-maximal) mineralization, and maximal mineralization at 14 days (18 days in total). Cells were harvested at 8 hours, 24 hours, 96 hours and 10 days for microarray profiling. Mineralization level was quantified at 14 days by Alizarin Red S staining (right column). GAGE was applied to infer the most differentially expressed pathways or gene sets between the matched samples with or without BMP6 at different time points. At time points corresponding to 8, 24 and 96 hours, GAGE compares between only two sample conditions. At the 10 day time point, GAGE compares between two mineralized conditions versus two non-mineralized conditions.

At data analysis level, we developed a novel temporal pathway analysis procedure. Gene set analysis (GSA) is a well established strategy to identify pathways or gene sets associated with a particular phenotype or condition [1316]. Temporary analysis was done on biological processes previously [17]. However, traditional GSA methods do not directly apply without externally specified selection criteria for temporal changes. Methods like GSEA [13] accept continuous phenotype labels hence can correlate gene set changes to time or pre-defined temporal patterns. To our knowledge, these methods do not infer temporal pathway level changes without reference patterns. For such temporal pathway analysis to be practical, a method needs to be: (1) applicable to datasets with small or even changing sample size at each time point or condition. (2) both sensitive and selective to capture subtle yet real regulatory signals over time. To meet these challenges, we designed a special analysis procedure (Figure 1) based on our newly developed GAGE (Generally Applicable Gene-set Enrichment) method [18].

Using this joint experimental and computational approach, we identified both the pathways and their temporal responses to BMP6 signaling during osteoblast differentiation and mineralization.


Following our previous study [5], we explored BMP6 induced human MSC osteoblast gene expression and function. Our preliminary experiments showed that 8 hours BMP6 exposure was sufficient to induce expression of osteoblast differentiation marker genes in human MSCs. A minimum of 24 hours exposure was required to form mineralized matrix at 14 days after the initiation of BMP treatment. A maximal mineralization response was observed upon 96 hours of BMP6 treatment.

We designed a microarray study to explore the regulatory mechanisms underneath these phenotypic changes at different stages of human MSC osteoblast differentiation (Figure 1). Our recently developed GAGE method [18] was applied to infer the significantly perturbed KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, GO (Gene Ontology) term groups, and experimentally derived gene sets (experimental sets for short) by BMP6 treatment at different times along the induction process (Figure 1). We examined these significant gene sets in details below.

Significantly perturbed KEGG pathways during BMP6 osteogenic induction

We observed a single set of KEGG pathways that were significantly perturbed at the gene expression level throughout the BMP6 induction process (Figure 2F and Table 1). This result suggests that these regulatory mechanisms are constantly involved at different stages of BMP6 induced osteoblast differentiation and mineralization. To test whether a KEGG pathway are significantly associated with a phenotype or a sample condition, we account for gene expression level perturbation in both directions (both up- and down-regulation), because a pathway commonly includes both positive and negative effectors and local feedback loops to keep the system in balance [18]. Therefore, we refer to a KEGG pathway as significantly perturbed rather than activated or inhibited as a whole.

Figure 2

The expression perturbation patterns induced by BMP6 treatment in eight significant KEGG pathways. These pathways are consistently significantly differentially expressed or near so based on GAGE. The mean t-statistics from multiple one-on-one comparisons between the two sample conditions is used as overall perturbation magnitude for each pathway. Perturbation magnitude here is defined as the absolute value of both positive and negative gene expression changes without considering the actual perturbation direction. The dashed line indicates a t-statistic of 1.79, which roughly corresponds to p = 0.01 for 8-96 h or 0.001 for 10 d. Panel (a) and (b) are divided only for visual clarity.

Table 1 Interpretation and validation of the significant KEGG pathways inferred by GAGE.

These significant KEGG pathways show different temporal perturbation patterns (Figure 2). The TGF-beta signaling pathway, Cytokine-cytokine receptor interaction, and Wnt signaling pathway are most perturbed at 8h. The Jak-STAT signaling pathway, MAPK signaling pathway and p53 signaling pathway are most perturbed at 24h. Focal adhesion and ECM-receptor interaction are most perturbed at 10d.

We explored the pattern of differential gene expression induced by BMP6 treatment in three representative pathways in detail (Figure 3 a-c and Additional file 1,2,3,4,5). First, TGF-beta signaling pathway is the top significant pathway (compared to other pathways) at 8h and is most perturbed (compared to itself at other time points) also at 8h. Therefore, this pathway is likely triggered directly by BMP6 treatment, and could be the signal initiating the osteoblast differentiation. This role of TGF-beta signaling is consistent with previous work [19] and the common sense: BMP6 as a canonical BMP triggers canonical BMP signaling, which is one major branch of TGF-beta signaling pathway. Second, focal adhesion is the most significant pathway (compared to other pathways) after 24h and is most perturbed (compared to itself at other time points) at 10 days. This pathway is likely the convergence point of the regulatory signals, and it is associated with late stage osteoblast differentiation and mineralization. These results are consistent with the role of focal adhesion inferred from experimental work in the literature [20, 21]. Third, the MAPK signaling pathway is the most perturbed at 24h or the middle stage of differentiation, which suggests that this pathway is likely an intermediate step during the BMP6 induced signal relay process at gene expression level. Indeed, MAPKs have been reported to mediate BMP effect during osteoblast differentiation [22, 23]. The temporal perturbation patterns in other significant pathways suggest their functions in the BMP6 induction process, which are supported by both observations in the literature and evidence from the expression data (Table 1). The temporal roles assigned to the pathways are relative, because all these pathways are important throughout BMP6 induced osteoblast differentiation and mineralization (Figure 2).

Figure 3

Gene expression fold changes induced by BMP6 in three representative significant KEGG pathways. Each pathway shown is at its most perturbed time point: (a) TGF-beta signaling pathway at 8h; (b) Focal adhesion at 10d; (c) MAPK signaling pathway at 24h. Gene expression level fold changes are standardized over the standard deviation of fold changes for all genes. The standardized fold changes are visualized by using KEGGanim web tool [56]. Note that one KEGG node may correspond to multiple closely related genes with the same function. Relevant other pathways are magnified locally for better view. Gene names are intentionally omitted by KEGGanim for clear view of the gene expression changes in pseudo-color.

Interestingly, the significant KEGG pathways are not distinct but rather act as an integrated super regulatory system. They interconnect to each other as shown on KEGG pathway graphs (Figure 3a-c). For example, the TGF-beta pathway triggers the MAPK signaling pathway (Figure 3a), whereas MAPK signaling pathway connects to focal adhesion (Figure 3b). These pathways also share common downstream response processes including apoptosis, cell cycle etc (Figure 3a-c). Using the top KEGG pathways inferred by GAGE [18], gene expression data and connections between pathways from LinkDB module of KEGG databases [24], we created an integrated dynamic network view of pathways in Figure 4. The upstream nodes including TGF-beta signaling pathway and cytokine-cytokine receptor interaction are most perturbed at the early stage. Downstream nodes including focal adhesion and ECM-receptor interaction are most perturbed at the late stage. Midstream nodes including MAPK signaling and p53 signaling are most perturbed at the middle stage. These sequential perturbation patterns across interconnected pathways suggest a dynamic transmission process of the regulatory signals induced by BMP6 treatment at the transcriptional level.

Figure 4

An integrated network of the significant KEGG pathways with their temporal perturbation patterns. Significant KEGG pathways (colored pie charts) or closely related other pathways (text only) are connected by arrows as indicated in the KEGG database. The mean t-statistics from multiple one-on-one comparisons at four different time points are plotted in pseudo heat color in the pie nodes as the average perturbation magnitude for significant KEGG pathways. Full names for some abbreviated KEGG pathways are: Cytokine-cytokine receptor interaction (Cytokine), Ubiquitin mediated proteolysis (Ublysis), Phosphatidylinositol signaling system (phosph).

The connected significant pathways (Figures 4) often share component genes that are evidently perturbed in expression (Table 2). Particularly, overlaps between significant KEGG pathways (Table 2) are consistent with the connections between pathways (Figure 4): high overlaps (Table 2) almost always suggest direct connections between pathways (Figure 4), and vise versa. Overlapping component genes serve as bridges across these relatively independent functional modules or pathways, hence perturbation in one pathway such as the BMP-TGF-beta signaling pathway can be propagated throughout other relevant pathways. These data suggest that these significant KEGG pathways work together as an integrated regulatory system.

Table 2 The overlaps in perturbed member genes between the significant KEGG pathways inferred by GAGE

Significantly perturbed GO term gene sets during BMP6 osteogenic induction

In contrast to the top KEGG pathways, the top significant GO term gene sets change with time (Figure 5 and Table 3). Most relevant GO term gene sets are significantly up or down regulated only for part of the induction process, except for the Notch signaling pathway and insulin-like growth factor receptor binding. In contrast to the KEGG pathways, where we test for two-directional perturbations, we test for one-directional changes in GO term gene sets as a complementary analysis to the KEGG pathway analysis above. More details are described in the Methods section. Additional file 6 includes comprehensive lists of significant KEGG pathways and GO groups at different time points.

Figure 5

The expression perturbation patterns induced by BMP6 treatment in nine significant GO term gene sets. Each of these GO term gene sets is significant in at least one time point based on GAGE. (a) GO:0007219 Notch signaling pathway, GO:0005520 insulin-like growth factor binding, GO:0005159 insulin-like growth factor receptor binding, GO:0008430 selenium binding, (b) GO:0007156 homophilic cell adhesion, GO:0006955 immune response, GO:0001501 skeletal development, GO:0001503 ossification, GO:0016491 oxidoreductase activity. The mean t-statistic from multiple one-on-one comparisons between the two sample conditions is used as a measure of the overall perturbation for each GO term. Different from KEGG pathways, perturbation direction is considered here. The dashed lines mark t = +/-1.30, which roughly correspond to p = 0.05 for 8-96 h or 0.01 for 10 d.

Table 3 Interpretation and validation information of the significant GO term gene sets inferred by GAGE

The top GO term gene sets suggest novel yet plausible regulatory processes involved in BMP6 induced osteoblast differentiation and mineralization. Some GO terms are self-explaining, including skeletal development and ossification (Figure 5b). Other less evident top GO terms are supported by the literature (Table 3). For example, the notch signaling pathway (Figure 5a) is directly involved in BMP2 induced osteoblast differentiation [11, 25], suggesting a possible connection to BMP6 too. Both insulin-like growth factor receptor binding and insulin-like growth factor binding (Figure 5a) suggest that the IGF signal is important for osteoblast differentiation and mineralization. Indeed IGF1 [26] and IGF1R [27] promote mineralization and bone formation, whereas IGFBP 3-6 [2831] sequesters IGF1 and inhibits osteoblast differentiation and mineralization. Direct connections between immune response (Figure 5b) and bone formation and metabolism have been well appreciated [32, 33], which merge into a new research area, osteoimmunology [33]. Genes from homophilic cell adhesion [34] and oxidoreductase activity [35] (Figure 5b) sets play a role in osteoblast differentiation. Selenium (Figure 5a) deficiency is associated with osteoporosis [36], a disease in bone formation and metabolism.

Significantly perturbed experimentally derived gene sets during BMP6 osteogenic induction

Our analysis also identified the most significantly up or down-regulated experimental sets during BMP6 induction. Experimental sets are gene groups coregulated or coexpressed in certain chemical or genetic perturbation experiments from the literature. A significant experiment set suggests that their common regulator was perturbed in our experiment too. There are commonly multiple experimental sets describing the same perturbation or mapping to the same regulatory mechanism. A non-redundant subset of the top experimental sets were collected in Table 4 with their perturbation patterns shown in Figure 6. Similar to gene sets based on GO terms but not to those on KEGG pathways, experimental sets were tested for perturbations in only one direction, either up- or down- regulated [18].

Table 4 Interpretation and validation information on the significant experimental sets inferred by GAGE
Figure 6

The expression perturbation patterns induced by BMP6 treatment in six significant experimental sets. Each of these experimental sets is significantly up or down regulated in at least one time point based on GAGE. (a) MYB targets, BAF57 down, BAF57 up; (b) IFNb up, VEGF up, IRS down. The mean t-statistics from multiple one-on-one comparisons between the two sample conditions is used as overall perturbation for each gene set. Different from KEGG pathways, perturbation direction is considered here. The dashed lines mark t = +/-1.78, which roughly correspond to p = 0.01 for 8-96 h or 0.001 for 10 d.

MYB transcription factor is identified as a novel regulator for osteoblast differentiation. MYB target gene set was down regulated at 8h, 24h and 10d with no change in MYB expression level. MYB transcriptional activity at protein level can be inhibited through two potential mechanisms following BMP6 treatment: the activation of Wnt signaling pathways (Table 1 and Figure 2) phosphorylates and degrades of MYB protein [37], and BMP/TGF-beta and Wnt responsive OVOL1 antagonizes transcriptional activation of MYB by competing for target promoter binding [38].

Another novel transcriptional regulator for osteoblast differentiation we predicted is BAF57, which is the regulatory subunit SWI/SNF chromatin remodeling complex [39]. Multiple BAF57 target genes are directly related to osteoblast differentiation and function (Additional file 1: Supplementary Table 1). Indeed, SWI/SNF regulates osteoblast-specific transcription through chromatin structure modification [40]. BMP6 treatment may target SWI/SNF to nucleus though SMAD1 signal [41] or p38 MAPK pathway targets SWI-SNF chromatin-remodeling complex [42]. Interestingly, BAF57 positive target genes are up-regulated and negative targets down-regulated during BMP6 induction, which further confirms the involvement of BAF57 activity. The different timing of the positive and negative regulation likely suggests different dynamics of these actions.

Other interesting regulatory mechanisms are inferred based on the top ranking experimental sets (Table 4). Interferon beta (interferon alpha and gamma too, but not shown) positive target gene sets are down-regulated, which is consistent with Jak-STAT pathway from KEGG (Table 1). VEGF positive targets are down-regulated, likely because VEGF gene expression is down-regulated due to MYB inhibition [43]. IRS negative targets are down-regulated is consistent with activation of IGF signal, particularly up-regulation of IGF receptor binding proteins (Table 4 and Figure 5).

Discussion and Conclusion

This is the first high throughput microarray study on BMP6 induced transcriptional program in human MSC. It covers the whole process from early to late stage osteoblast differentiation and mineralization. We conducted a comprehensive gene set analysis to identify relevant regulatory mechanisms and functional groups. We inferred a series of significant KEGG pathways, GO terms and experimental sets at different stages of BMP6 induction process. We not only showed which pathways or gene sets are significant, but also when and how they are involved in the osteoblast differentiation and mineralization. Different from common pathway analyses [13, 14, 16], our work further captures the interconnections among individual pathways or functional groups and integrate them into a whole system. Taken together, this work provides clearer mechanistic picture of osteoblast differentiation and function.

We inferred novel and coherent sets of regulatory mechanisms downstream of BMP6 signaling during osteoblast differentiation and mineralization. First, the same set of KEGG pathways are constantly involved in BMP6 induction. Their roles in osteogenic induction are clarified based on their perturbation patterns and connected to relevant discoveries in literature. These significant KEGG pathways are not separated but rather they work as a unified super regulatory system, and the pathway perturbation patterns we derived reflect a dynamic transmission process of the regulatory signal at transcriptional level along the super system. Second, a varying set of GO processes and functional groups are involved at different stage of BMP6 induced osteoblast differentiation and mineralization. These suggest novel yet plausible regulatory mechanisms, which are connected to but have not directly and explicitly introduced in literature works. Third, the most significant experimental sets suggest novel transcriptional regulators including MYB and BAF57, and regulatory pathways consistent with predictions based on KEGG and GO gene sets above.

Connections between KEGG pathways are evident as shown in the super regulatory network of pathways (Figure 4 and Table 2). Perturbations propagate along the super network at two levels: at protein level, the phosphorylation, binding, activation/inhibition events relay along pathways and transmit into interconnected pathways, as stated by KEGG graphs (Figure 3); at transcriptional level, gene expression perturbation propagates through auto-regulatory (feedback and feed forward) loops within pathways, and bridges into its neighbor pathways through the multiple component genes in common, as suggested by our pathway analysis results (Figure 4 and Table 2). Protein level transmission is faster, but transcriptional level transmission lasts longer hence ensure the long term biological effects. These longer lasting transcriptional effects are clearly evident in our data as the gene expression levels were perturbed long after the withdrawal of BMP6 treatment. The hard wired KEGG pathways and interconnections between them define how BMP6 signal triggers downstream programs (Figure 4).

There are also connections between BMP6 signal and the processes/groups represented by GO terms and experimental sets. For example, Notch signal and IGF signal are involved in the whole induction process (Figure 5a) like all significant KEGG pathways (Figure 2). It follows that these two signals should be also part of the super regulatory system and interconnected with multiple significant KEGG pathways of the network (Figure 4). This hypothesis is well supported by our data and literature: (1) Notch signal directly interacts with BMP signal. SMAD1 and NIC synergize to induce expression of HEY1 and other Notch targets [4446]. Indeed up-regulation of HEY1 (and HEY2) requires continued BMP6 treatment by 96 hours (Additional file 1: Supplementary Figure 4c). Besides this binding-synergy at protein level, Notch ligand JAG1 expression is also up-regulated directly by BMP6 (Additional file 1: Supplementary Figure 4c). Notch also interacts with Wnt signal [47]. (2) Growth hormone (GH) signal [48] as part of cytokine-cytokine receptor interaction (KEGG pathway 04060) and Jak-STAT pathways (KEGG pathway 04630) are activated and by BMP signal (Figure 4) with GHR up-regulated (Additional file 1: Supplementary Figure 4a-b), GH signal up-regulates (IGF1 and IRS1, Additional file 1: Supplementary Figure 4d) and activates IGF signal in turn [48]. Connections between BMP signal and the two predicted transcriptional regulators, MYB and BAF57, are described above in the Results.

When examined together, we find a consistent picture emerging from the lists of significant KEGG pathways, GO terms, and experimental sets. For example, KEGG focal adhesion and ECM receptor interaction pathways (Table 1), GO homophilic cell adhesion (Table 3) and extracellular matrix structural constituent (not shown) groups consistently show the relevance of cell adhesion and extracellular matrix in osteoblast differentiation and mineralization. GO immune response groups (Table 3) echoes KEGG cytokine-cytokine receptor interaction pathway (Table 1). Similarly, significant experimental target genes sets (Table 4) closely reflected changes in the regulatory KEGG pathways (Table 1) or GO processes/groups (Table 3).

Interestingly, there are discrepancies among the clustering and the significant KEGG pathways, GO terms, and experimental sets. For example, Notch signaling is defined as both a KEGG pathway and a GO process. This KEGG pathway is not significant (not shown) but this GO process is (Figure 5). This discrepancy arises from two sources: (1) different definitions, i.e. KEGG pathways contain partially different set of genes from corresponding GO processes. While KEGG pathways tend to cover the whole homeostatic signal transmission systems even across multiple transcriptional cycles, GO usually covers one or multiple discrete steps or functional groups for a process. KEGG and GO definitions can be considered complementary and both provide valuable gene sets for our analysis. (2) GAGE [18] treats KEGG pathways and GO term gene sets differently: genes under a GO term are taken as a group coregulated towards a single direction, either all up or all down regulated, whereas genes in a KEGG pathway are frequently not coregulated and expression changes in both directions are counted. Timing discrepancies exist between experimental sets and corresponding KEGG pathways. For example, IFN positive target sets (only IFN beta shown) are not significant at 24-96 hours and MYB target set not at 96 hours (Figure 6) while Jak-STAT pathway and Wnt signaling pathway are significantly perturbed all the time (Figure 2). This can be explained by the fact that two-directional perturbation treatment for KEGG pathway does not account for direction or net effect of the perturbation, whether inhibited, activated or no overall effect. In the other hand, GO term analysis has no such issue, and IRS negative set and corresponding GO IGF receptor binding group are both significant all the time.

In this work, we took a systems approach in studying MSC differentiation. We combined experimental and computational work to reconstruct a unified picture of BMP6 signaling. The same set of experimental design and computational approach could be used to study other physiological and pathological processes, such as the differentiation of other cell types or tumorogenesis.


Cell culture and BMP6 osteogenic induction

Passage 5 human MSC (5 × 105) were plated in 24-well dishes and cultured for 3 days. The cells were subsequently placed in serum free media supplemented with ITS for 24 hours. BMP6 was then added for the pre-defined time periods and then removed as shown in Figure 1. To remove BMP6, cells were rinsed 2 times before fresh media without BMP6 was added. Ascorbate and b-glycerolphosphate were added 4 days after the initiation of BMP treatment. Cells were harvested at the indicated time points. To quantifying mineralization, the plates were stained with Alizarin Red S [5] 14 days after the initiation of BMP treatment.

Microarray experiment and analysis

Human MSC cells underwent osteogenic induction with BMP6 treatment for 0 hours, 8 hours, 24 hours, and 96 hours, which correspond to four phenotypic groups, i.e. control, preosteoblast (no mineralization), (sub-maximal) mineralization, and maximal mineralization at 14 days after the initiation of BMP treatment (18 days in total). Cells were harvested at 8 hours, 24 hours, 96 hours and 10 days for microarray profiling using Affymetrix U133 plus GeneChip® platform. Assays were run in duplicate for a total of 20 arrays.

The raw data were processed by using FRAMS [49] with up-to-date probe set definition [44]. GAGE [18] was applied to infer the most differentially expressed pathways or gene sets between the BMP6 added or withdrawn samples under comparison at different time points (Figure 1). As is shown in Figure 1, the comparisons at 8, 24 and 96 hours were between two sample conditions, whereas comparison at 10 days was between the 24 and 96 hours BMP6 groups versus 0 and 8 hours BMP6 groups.

Pathway analysis using GAGE

We use GAGE [18], Generally Applicable Gene set Enrichment, a novel method we developed for gene set/pathway analysis. The GAGE method is implemented in R in the "gage" package, available through Bioconductor at As the GAGE procedure has been described in detail in the original paper [18], here is a brief summary of the method.

Step 1: Gene sets separation. Gene sets are derived or collected from KEGG pathways [24], GO [50] and MSigDB [51] databases, as KEGG pathways, GO terms and experimental sets respectively. GAGE treats KEGG pathways differently from GO terms and experimental sets: member genes for a GO term or experimental set are taken as a group coregulated towards a single direction, either all up or all down regulated, whereas genes in a KEGG pathway are frequently not coregulated and expression changes in both directions are counted. We test for one-directional changes in GO terms because common function is frequently connected to coregulation and coexpression [52, 53], which has also been extensively used in classical GO analysis tools [13, 54, 55]. Such treatment will miss some pathway-like GO terms which are significantly perturbed in two directions, but we intended to have GO analysis as a complementary work to the KEGG pathway analysis. Nonetheless, GAGE provides two-directional test options for all types of gene sets. For GO analysis, all terms are included without differentiating categories, i.e. Biological Process, Cellular Component, Molecular Function, and without considering the hierarchical organization of the ontology tree.

Step 2: One-on-one comparisons. Instead of comparing BMP6 treated samples vs controls as two groups, GAGE does one-on-one comparison between samples from the two groups at a time. For each one-on-one comparison, log based fold changes are calculated for all genes. GAGE conducts two-sample t-test on the average fold change in specific gene sets against that for the background of the whole set. This one-on-one comparison procedure is repeated for all potential experiment-control pairs.

Step 3: Summarization. For each gene set, GAGE derives a global P-value based on a meta test on the negative log sum of multiple P-values for this set from all one-on-one comparisons between experiments and controls.

Overlap between significant pathways

Significant KEGG pathways may overlap in terms of perturbed member genes, which could suggest biologically meaningful connections between pathways. A gene is counted as perturbed when its absolute fold change is at least one standard deviation higher than the mean of all genes. The significance of overlap between pathways is inferred using hypergeometric test as following. Let n represent number of perturbed genes in a pathway, N number of all perturbed genes in a two-state comparison, x the overlap perturbed genes between pathways, pathways 1 and 2 are labeled by corresponding suffix. Then the chance to see X or more overlap perturbed genes between two significant pathways is:


Perturbation pattern visualization

For significant KEGG pathways, we generated graphs to visualize the dynamic expression perturbation at two levels: individual genes (Figure 3 and Additional file 2,3,4,5) and whole pathways (Figure 4). Gene expression level fold changes are standardized over the standard deviation of fold changes for all genes. The standardized fold changes for individual genes in KEGG pathways are visualized by using KEGGanim web tool [56] in Figure 3 and Additional file 2,3,4,5. We present an integrated network to show the connections between pathways and average expression perturbations for them in Figure 4. Significant KEGG pathways or closely related other pathways connected to them are represented by nodes. Connections between these pathways collected from KEGG database and graphs [24] are represented by arrows plus edges. The mean t-statistics from two-sample t-test from multiple one-on-one comparisons are plotted in pseudo heat colors as average perturbation magnitude for significant KEGG pathways.


  1. 1.

    Harada S, Rodan GA: Control of osteoblast function and regulation of bone mass. Nature. 2003, 423: 349-55. 10.1038/nature01660

  2. 2.

    Chen D, Zhao M, Mundy GR: Bone morphogenetic proteins. Growth Factors. 2004, 22: 233-41. 10.1080/08977190412331279890

  3. 3.

    Li X, Cao X: BMP signaling and skeletogenesis. Ann N Y Acad Sci. 2006, 1068: 26-40. 10.1196/annals.1346.006

  4. 4.

    Lavery KS, Swain PM, Falb D, Alaoui-Ismaili MH: BMP-2/4 and BMP-6/7 differentially utilize cell surface receptors to induce osteoblastic differentiation of human bone marrow derived mesenchymal stem cells. J Biol Chem. 2008

  5. 5.

    Friedman MS, Long MW, Hankenson KD: Osteogenic differentiation of human mesenchymal stem cells is regulated by bone morphogenetic protein-6. J Cell Biochem. 2006, 98: 538-54. 10.1002/jcb.20719

  6. 6.

    de Jong DS, Vaes BL, Dechering KJ, Feijen A, Hendriks JM, Wehrens R, Mummery CL, van Zoelen EJ, Olijve W, Steegenga WT: Identification of novel regulators associated with early-phase osteoblast differentiation. J Bone Miner Res. 2004, 19: 947-58. 10.1359/JBMR.040216

  7. 7.

    Peng Y, Kang Q, Cheng H, Li X, Sun MH, Jiang W, Luu HH, Park JY, Haydon RC, He TC: Transcriptional characterization of bone morphogenetic proteins (BMPs)-mediated osteogenic signaling. J Cell Biochem. 2003, 90: 1149-65. 10.1002/jcb.10744

  8. 8.

    Harris SE, Guo D, Harris MA, Krishnaswamy A, Lichtler A: Transcriptional regulation of BMP-2 activated genes in osteoblasts using gene expression microarray analysis: role of Dlx2 and Dlx5 transcription factors. Front Biosci. 2003, 8: s1249-65. 10.2741/1170

  9. 9.

    Korchynskyi O, Dechering KJ, Sijbers AM, Olijve W, ten Dijke P: Gene array analysis of bone morphogenetic protein type I receptor-induced osteoblast differentiation. J Bone Miner Res. 2003, 18: 1177-85. 10.1359/jbmr.2003.18.7.1177

  10. 10.

    Kalajzic I, Staal A, Yang WP, Wu Y, Johnson SE, Feyen JH, Krueger W, Maye P, Yu F, Zhao Y, Kuo L, Gupta RR, Achenie LE, Wang HW, Shin DG, Rowe DW: Expression profile of osteoblast lineage at defined stages of differentiation. J Biol Chem. 2005, 280: 24618-26. 10.1074/jbc.M413834200

  11. 11.

    de Jong DS, Steegenga WT, Hendriks JM, van Zoelen EJ, Olijve W, Dechering KJ: Regulation of Notch signaling genes during BMP2-induced differentiation of osteoblast precursor cells. Biochem Biophys Res Commun. 2004, 320: 100-7. 10.1016/j.bbrc.2004.05.150

  12. 12.

    Balint E, Lapointe D, Drissi H, van der Meijden C, Young DW, van Wijnen AJ, Stein JL, Stein GS, Lian JB: Phenotype discovery by gene expression profiling: mapping of biological processes linked to BMP-2-mediated osteoblast differentiation. J Cell Biochem. 2003, 89: 401-26. 10.1002/jcb.10515

  13. 13.

    Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005, 102: 15545-50. 10.1073/pnas.0506580102

  14. 14.

    Kelder T, Conklin BR, Evelo CT, Pico AR: Finding the right questions: exploratory pathway analysis to enhance biological discovery in large datasets. PLoS Biol. 2010, 8: e1000472- 10.1371/journal.pbio.1000472

  15. 15.

    Davies MN, Meaburn EL, Schalkwyk LC: Gene set enrichment; a problem of pathways. Brief Funct Genomics. 9: 385-90.

  16. 16.

    Nam D, Kim SY: Gene-set approach for expression pattern analysis. Brief Bioinform. 2008, 9: 189-97. 10.1093/bib/bbn001

  17. 17.

    Kleemann R, van Erk M, Verschuren L, van den Hoek AM, Koek M, Wielinga PY, Jie A, Pellis L, Bobeldijk-Pastorova I, Kelder T, Toet K, Wopereis S, Cnubben N, Evelo C, van Ommen B, Kooistra T: Time-resolved and tissue-specific systems analysis of the pathogenesis of insulin resistance. PLoS One. 5: e8817-

  18. 18.

    Luo W, Friedman MS, Shedden K, Hankenson KD, Woolf PJ: GAGE: generally applicable gene set enrichment for pathway analysis. BMC Bioinformatics. 2009, 10: 161- 10.1186/1471-2105-10-161

  19. 19.

    Gosselet FP, Magnaldo T, Culerrier RM, Sarasin A, Ehrhart JC: BMP2 and BMP6 control p57(Kip2) expression and cell growth arrest/terminal differentiation in normal primary human epidermal keratinocytes. Cell Signal. 2007, 19: 731-9. 10.1016/j.cellsig.2006.09.006

  20. 20.

    Kim JB, Leucht P, Luppen CA, Park YJ, Beggs HE, Damsky CH, Helms JA: Reconciling the roles of FAK in osteoblast differentiation, osteoclast remodeling, and bone regeneration. Bone. 2007, 41: 39-51. 10.1016/j.bone.2007.01.024

  21. 21.

    Salasznyk RM, Klees RF, Williams WA, Boskey A, Plopper GE: Focal adhesion kinase signaling pathways regulate the osteogenic differentiation of human mesenchymal stem cells. Exp Cell Res. 2007, 313: 22-37. 10.1016/j.yexcr.2006.09.013

  22. 22.

    Ulsamer A, Ortuno MJ, Ruiz S, Susperregui AR, Osses N, Rosa JL, Ventura F: BMP-2 induces Osterix expression through up-regulation of Dlx5 and its phosphorylation by p38. J Biol Chem. 2008, 283: 3816-26.

  23. 23.

    Celil AB, Campbell PG: BMP-2 and insulin-like growth factor-I mediate Osterix (Osx) expression in human mesenchymal stem cells via the MAPK and protein kinase D signaling pathways. J Biol Chem. 2005, 280: 31353-9. 10.1074/jbc.M503845200

  24. 24.

    KEGG: Kyoto Encyclopedia of Genes and Genomes.

  25. 25.

    Nobta M, Tsukazaki T, Shibata Y, Xin C, Moriishi T, Sakano S, Shindo H, Yamaguchi A: Critical regulation of bone morphogenetic protein-induced osteoblastic differentiation by Delta1/Jagged1-activated Notch1 signaling. J Biol Chem. 2005, 280: 1584-8.

  26. 26.

    Zhao G, Monier-Faugere MC, Langub MC, Geng Z, Nakayama T, Pike JW, Chernausek SD, Rosen CJ, Donahue LR, Malluche HH, Fagin JA, Clemens TL: Targeted overexpression of insulin-like growth factor I to osteoblasts of transgenic mice: increased trabecular bone volume without increased osteoblast proliferation. Endocrinology. 2000, 141: 2674-82. 10.1210/en.141.7.2674

  27. 27.

    Zhang M, Xuan S, Bouxsein ML, von Stechow D, Akeno N, Faugere MC, Malluche H, Zhao G, Rosen CJ, Efstratiadis A, Clemens TL: Osteoblast-specific knockout of the insulin-like growth factor (IGF) receptor gene reveals an essential role of IGF signaling in bone matrix mineralization. J Biol Chem. 2002, 277: 44005-12. 10.1074/jbc.M208265200

  28. 28.

    Strohbach C, Kleinman S, Linkhart T, Amaar Y, Chen ST, Mohan S, Strong D: Potential involvement of the interaction between insulin-like growth factor binding protein (IGFBP)-6 and LIM mineralization protein (LMP)-1 in regulating osteoblast differentiation. J Cell Biochem. 2008

  29. 29.

    Mukherjee A, Rotwein P: Insulin-like growth factor-binding protein-5 inhibits osteoblast differentiation and skeletal growth by blocking insulin-like growth factor actions. Mol Endocrinol. 2008, 22: 1238-50. 10.1210/me.2008-0001

  30. 30.

    Silha JV, Mishra S, Rosen CJ, Beamer WG, Turner RT, Powell DR, Murphy LJ: Perturbations in bone formation and resorption in insulin-like growth factor binding protein-3 transgenic mice. J Bone Miner Res. 2003, 18: 1834-41. 10.1359/jbmr.2003.18.10.1834

  31. 31.

    Zhang M, Faugere MC, Malluche H, Rosen CJ, Chernausek SD, Clemens TL: Paracrine overexpression of IGFBP-4 in osteoblasts of transgenic mice decreases bone turnover and causes global growth retardation. J Bone Miner Res. 2003, 18: 836-43. 10.1359/jbmr.2003.18.5.836

  32. 32.

    Arron JR, Choi Y: Bone versus immune system. Nature. 2000, 408: 535-6. 10.1038/35046196

  33. 33.

    Walsh MC, Kim N, Kadono Y, Rho J, Lee SY, Lorenzo J, Choi Y: Osteoimmunology: interplay between the immune system and bone metabolism. Annu Rev Immunol. 2006, 24: 33-63. 10.1146/annurev.immunol.24.021605.090646

  34. 34.

    Arai F, Ohneda O, Miyamoto T, Zhang XQ, Suda T: Mesenchymal stem cells in perichondrium express activated leukocyte cell adhesion molecule and participate in bone marrow formation. J Exp Med. 2002, 195: 1549-63. 10.1084/jem.20011700

  35. 35.

    Bennett KP, Bergeron C, Acar E, Klees RF, Vandenberg SL, Yener B, Plopper GE: Proteomics reveals multiple routes to the osteogenic phenotype in mesenchymal stem cells. BMC Genomics. 2007, 8: 380- 10.1186/1471-2164-8-380

  36. 36.

    Ebert R, Jakob F: Selenium deficiency as a putative risk factor for osteoporosis. International Congress Series: Nutritional Aspects of Osteoporosis 2006. Proceedings of the 6th International Symposium on Nutritional Aspects of Osteoporosis, 4-6 May 2006, Lausanne, Switzerland. 2007, 1297: 158-164.

  37. 37.

    Kanei-Ishii C, Ninomiya-Tsuji J, Tanikawa J, Nomura T, Ishitani T, Kishida S, Kokura K, Kurahashi T, Ichikawa-Iwata E, Kim Y, Matsumoto K, Ishii S: Wnt-1 signal induces phosphorylation and degradation of c-Myb protein via TAK1, HIPK2, and NLK. Genes Dev. 2004, 18: 816-29. 10.1101/gad.1170604

  38. 38.

    Nair M, Bilanchone V, Ortt K, Sinha S, Dai X: Ovol1 represses its own transcription by competing with transcription activator c-Myb and by recruiting histone deacetylase activity. Nucleic Acids Res. 2007, 35: 1687-97. 10.1093/nar/gkl1141

  39. 39.

    Garcia-Pedrero JM, Kiskinis E, Parker MG, Belandia B: The SWI/SNF chromatin remodeling subunit BAF57 is a critical regulator of estrogen receptor function in breast cancer cells. J Biol Chem. 2006, 281: 22656-64. 10.1074/jbc.M602561200

  40. 40.

    Villagra A, Cruzat F, Carvallo L, Paredes R, Olate J, van Wijnen AJ, Stein GS, Lian JB, Stein JL, Imbalzano AN, Montecino M: Chromatin remodeling and transcriptional activity of the bone-specific osteocalcin gene require CCAAT/enhancer-binding protein beta-dependent recruitment of SWI/SNF activity. J Biol Chem. 2006, 281: 22695-706. 10.1074/jbc.M511640200

  41. 41.

    Colland F, Jacq X, Trouplin V, Mougin C, Groizeleau C, Hamburger A, Meil A, Wojcik J, Legrain P, Gauthier JM: Functional proteomics mapping of a human signaling pathway. Genome Res. 2004, 14: 1324-32. 10.1101/gr.2334104

  42. 42.

    Simone C, Forcales SV, Hill DA, Imbalzano AN, Latella L, Puri PL: p38 pathway targets SWI-SNF chromatin-remodeling complex to muscle-specific loci. Nat Genet. 2004, 36: 738-43. 10.1038/ng1378

  43. 43.

    Lutwyche JK, Keough RA, Hunter J, Coles LS, Gonda TJ: DNA binding-independent transcriptional activation of the vascular endothelial growth factor gene (VEGF) by the Myb oncoprotein. Biochem Biophys Res Commun. 2006, 344: 1300-7. 10.1016/j.bbrc.2006.04.045

  44. 44.

    Itoh F, Itoh S, Goumans MJ, Valdimarsdottir G, Iso T, Dotto GP, Hamamori Y, Kedes L, Kato M, ten Dijke Pt P: Synergy and antagonism between Notch and BMP receptor signaling pathways in endothelial cells. Embo J. 2004, 23: 541-51. 10.1038/sj.emboj.7600065

  45. 45.

    Dahlqvist C, Blokzijl A, Chapman G, Falk A, Dannaeus K, Ibanez CF, Lendahl U: Functional Notch signaling is required for BMP4-induced inhibition of myogenic differentiation. Development. 2003, 130: 6089-99. 10.1242/dev.00834

  46. 46.

    Takizawa T, Ochiai W, Nakashima K, Taga T: Enhanced gene activation by Notch and BMP signaling cross-talk. Nucleic Acids Res. 2003, 31: 5723-31. 10.1093/nar/gkg778

  47. 47.

    Deregowski V, Gazzerro E, Priest L, Rydziel S, Canalis E: Notch 1 overexpression inhibits osteoblastogenesis by suppressing Wnt/beta-catenin but not bone morphogenetic protein signaling. J Biol Chem. 2006, 281: 6203-10. 10.1074/jbc.M508370200

  48. 48.

    Leroith D, Nissley P: Knock your SOCS off!. J Clin Invest. 2005, 115: 233-6.

  49. 49.

    Hochreiter S, Clevert DA, Obermayer K: A new summarization method for Affymetrix probe level data. Bioinformatics. 2006, 22: 943-9. 10.1093/bioinformatics/btl033

  50. 50.

    Gene Ontology.

  51. 51.

    The Molecular Signature Database.

  52. 52.

    Horan K, Jang C, Bailey-Serres J, Mittler R, Shelton C, Harper JF, Zhu JK, Cushman JC, Gollery M, Girke T: Annotating genes of known and unknown function by large-scale coexpression analysis. Plant Physiol. 2008, 147: 41-57. 10.1104/pp.108.117366

  53. 53.

    Allocco DJ, Kohane IS, Butte AJ: Quantifying the relationship between co-expression, co-regulation and gene function. BMC Bioinformatics. 2004, 5: 18- 10.1186/1471-2105-5-18

  54. 54.

    Al-Shahrour F, Diaz-Uriarte R, Dopazo J: FatiGO: a web tool for finding significant associations of Gene Ontology terms with groups of genes. Bioinformatics. 2004, 20: 578-80. 10.1093/bioinformatics/btg455

  55. 55.

    Zhang B, Schmoyer D, Kirov S, Snoddy J: GOTree Machine (GOTM): a web-based platform for interpreting sets of interesting genes using Gene Ontology hierarchies. BMC Bioinformatics. 2004, 5: 16- 10.1186/1471-2105-5-16

  56. 56.

    Adler P, Reimand J, Janes J, Kolde R, Peterson H, Vilo J: KEGGanim: pathway animations for high-throughput data. Bioinformatics. 2008, 24: 588-90. 10.1093/bioinformatics/btm581

  57. 57.

    Hughes FJ, Collyer J, Stanfield M, Goodman SA: The effects of bone morphogenetic protein-2, -4, and -6 on differentiation of rat osteoblast cells in vitro. Endocrinology. 1995, 136: 2671-7. 10.1210/en.136.6.2671

  58. 58.

    Takayanagi H, Sato K, Takaoka A, Taniguchi T: Interplay between interferon and other cytokine systems in bone metabolism. Immunol Rev. 2005, 208: 181-93. 10.1111/j.0105-2896.2005.00337.x

  59. 59.

    Fujita K, Janz S: Attenuation of WNT signaling by DKK-1 and -2 regulates BMP2-induced osteoblast differentiation and expression of OPG, RANKL and M-CSF. Mol Cancer. 2007, 6: 71- 10.1186/1476-4598-6-71

  60. 60.

    Rawadi G, Vayssiere B, Dunn F, Baron R, Roman-Roman S: BMP-2 controls alkaline phosphatase expression and osteoblast mineralization by a Wnt autocrine loop. J Bone Miner Res. 2003, 18: 1842-53. 10.1359/jbmr.2003.18.10.1842

  61. 61.

    Rajan P, Panchision DM, Newell LF, McKay RD: BMPs signal alternately through a SMAD or FRAP-STAT pathway to regulate fate choice in CNS stem cells. J Cell Biol. 2003, 161: 911-21. 10.1083/jcb.200211021

  62. 62.

    Wang X, Kua HY, Hu Y, Guo K, Zeng Q, Wu Q, Ng HH, Karsenty G, de Crombrugghe B, Yeh J, Li B: p53 functions as a negative regulator of osteoblastogenesis, osteoblast-dependent osteoclastogenesis, and bone remodeling. J Cell Biol. 2006, 172: 115-25. 10.1083/jcb.200507106

  63. 63.

    Lengner CJ, Steinman HA, Gagnon J, Smith TW, Henderson JE, Kream BE, Stein GS, Lian JB, Jones SN: Osteoblast differentiation and skeletal development are regulated by Mdm2-p53 signaling. J Cell Biol. 2006, 172: 909-21. 10.1083/jcb.200508130

  64. 64.

    Xiao G, Gopalakrishnan R, Jiang D, Reith E, Benson MD, Franceschi RT: Bone morphogenetic proteins, extracellular matrix, and mitogen-activated protein kinase signaling pathways are required for osteoblast-specific gene expression and differentiation in MC3T3-E1 cells. J Bone Miner Res. 2002, 17: 101-10. 10.1359/jbmr.2002.17.1.101

  65. 65.

    Xiao G, Wang D, Benson MD, Karsenty G, Franceschi RT: Role of the alpha2-integrin in osteoblast-specific gene expression and activation of the Osf2 transcription factor. J Biol Chem. 1998, 273: 32988-94. 10.1074/jbc.273.49.32988

Download references


This work was supported by NIH grant R01 DE017471. In addition, PJW is supported by NIH grant U54-DA-021519 and KDH is supported by NIH grants R01 AR054714 and R01 AR049682.

Author information



Corresponding authors

Correspondence to Kurt D Hankenson or Peter J Woolf.

Additional information

Authors' contributions

WL, MSF, KDH and PJW conceived and designed the study; WL conducted the computational analysis and result interpretation; MSF conducted the microarray experiment. WL and PJW drafted the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

Authors’ original submitted files for images

Rights and permissions

Reprints and Permissions

About this article

Cite this article

Luo, W., Friedman, M.S., Hankenson, K.D. et al. Time series gene expression profiling and temporal regulatory pathway analysis of BMP6 induced osteoblast differentiation and mineralization. BMC Syst Biol 5, 82 (2011).

Download citation


  • Bone Morphogenetic Protein
  • Osteoblast Differentiation
  • KEGG Pathway
  • Bone Morphogenetic Protein Signal
  • BMP6 Treatment