Similarity-based metric analysis approach for predicting osteogenic differentiation correlation coefficients and discovering the novel osteogenic-related gene FOXA1 in BMSCs

Data collection

The human PPI network including 19,344 genes is downloaded from the STRING database (Szklarczyk et al., 2016). It is then simplified by removing multiple edges and self-loops. Osteogenic differentiation-related genes are collected from online databases including GO, KEGG, BIOCARTA, PID, REACTOME, and WikiPathways. Among these databases, osteogenic differentiation-related gene sets or pathways are first retrieved using the keyword osteoblast. As a result, 10 osteogenic differentiation-related pathways including 337 unique genes are collected in our analysis. To determine the interaction between FOXA1 and osteogenic differentiation, four signaling pathways including hedgehog signaling pathway, BMP signaling pathway, ERK pathway, and Wnt signaling that have widely been studied to be associated with osteogenic differentiation are retrieved from The Molecular Signatures Databases (MSigDB) (Subramanian et al., 2005).

Construction of similarity matrix between different genes

To calculate the similarity matrix between genes, we proposed two indexes Jaccard similarity and Sorensen-Dice similarity to measure the topological relevance of genes in the human PPI network. Given a graph $G\left(V,E\right)$, Jaccard similarity between any two vertexes can be defined as the following formula: (1)${\displaystyle J\left(A,B\right)=\frac{|N\left(A\right)\cap N\left(B\right)|}{|N\left(A\right)\cup N\left(B\right)|}}$ where N(A) is the neighbors of vertex A and N(B) is the neighbors of vertex B. It calculates the proportion of co-neighbors among all neighbors of vertexes A and B in which a big similarity score indicates significant overlap of neighbors. In that way, the topological similarity between any two genes can be caught and extracted by evaluating the overlap of co-neighbors.

Additionally, another similarity-based index Sorensen-Dice similarity is carried out to measure the topological similarity of genes which can be defined as: (2)${\displaystyle SI\left(A,B\right)=\frac{2|N\left(A\right)\cap N\left(B\right)|}{\left|N\left(A\right)\right|+\left|N\left(B\right)\right|}.}$

Similar with Jaccard similarity, Sorensen-Dice similarity can utilize and calculate the topological similarity of genes by evaluating the topological overlap of neighbors. But it is normalized by dividing the total neighbors.

After evaluating the topological similarity of genes in the human PPI network based on Jaccard similarity and Sorensen-Dice similarity, these two indexes are then integrated them into a final similarity matrix which is defined as: (3)${\displaystyle S\left(A,B\right)=\frac{J\left(A,B\right)+SI\left(A,B\right)}{2}.}$

It takes into consideration both Jaccard similarity and Sorensen-Dice similarity and calculates the average similarity which indicates the affinity of genes. This similarity matrix is then used for osteo-specific scoring in the following analysis. All of these analyses are conducted using the igraph R package.

Functional enrichment analysis

To check whether these candidate genes identified by similarity index are enriched in some specific biological functions and states, we performed functional enrichment analysis of these genes in the GO and KEGG databases separately. We set 0.01 as the statistically significant threshold for BH-adjusted p-value. Pathways with BH-adjusted p-value less than 0.01 are considered to be statistically significant and presented by these candidate genes. GO enrichment analysis is performed in three different categories including biological process, cellular component, and molecular function. All of these analyses including functional enrichment analysis and visualization are implemented using clusterProfiler and enrichplot R packages (Yu et al., 2012).

Permutation test

To check whether these osteogenic differentiation-related pathways can be regulated by FOXA1, we performed a permutation test to investigate the similarity score with four well-known osteogenic differentiation-related pathways including hedgehog signaling pathway, BMP signaling, ERK pathway, and Wnt signaling pathway. As a reference, we select 50 genes randomly from the whole gene list to calculate the average similarity score with FOXA1. This sampling process is repeated 10,000 times. Furthermore, the interactions of FOXA1 with these four pathways are also evaluated by averaging the similarity score with these genes involved in the osteogenic pathway. The similarity score between osteogenic differentiation-related pathways and those randomly selected gene sets are then compared to determine whether FOXA1 is significantly associated with these pathways which indicates a potential regulation mechanism during osteogenic differentiation.

Cell culture and reagents

Cyagen Biosciences supplied hBMSCs (human bone mesenchymal stem cells, product number: HUXMA-01001, Cyagen Biosciences, Guangzhou, China), which may develop into osteoblasts, chondrocytes, and adipocytes under certain inductive circumstances. Adherent hBMSCs were cultured in culture flasks in a specific complete growth medium (HUXMA-90011; Cyagen Biosciences, Inc., Guangzhou, China) in a cell incubator at 37 ° C with 5% CO2 and passaged at around 80–90% confluence. Cells from passages two through six were utilized in future investigations.

Lentiviral packaging and cell infection

Obio Technology (Shanghai, China) provided the hBMSCs with a lentiviral package that included lentiviral particles to overexpress (FOXA1 overexpress group, OE) and overexpress control particles (FOXA1 overexpress negative control group, OE-NC), knock down FOXA1 (FOXA1 knockdown group, KD), and knockdown control particles. When hBMSCs achieved 30–50% confluence, lentiviral particles containing 5 ug/ml polybrene were introduced to the growing medium per the manufacturer’s instructions.

ALP staining and quantitation

Cells grown in osteogenic induction media for 5 days in 24-well plates were washed three times with phosphate-buffered saline, then fixed with 4% paraformaldehyde (BOSTER, Wuhan, China) for 30 min at room temperature. The cells were stained with a BCIP/NBT ALP colour development kit (Beyotime). To evaluate ALP activity, an ALP activity assay (BOSTER, Wuhan, China) was used in accordance with product instructions.

Alizarin Red S staining and quantitation

HBMSCs were cultivated in osteogenic induction media for 16 days after being passaged on 24-well plates. The cells were fixed in 4% paraformaldehyde (BOSTER, Wuhan, China) for 30 min at room temperature after being washed three times with PBS. Alizarin Red S solution (Cyagen Biosciences, Guangzhou, China)) was then added, and incubated at room temperature for 15 min. ARS stain was incubated with 10% cetylpyridinium chloride for 1 h at room temperature, then the solutions were collected, plated on a 96-well plate and measured at 560 nm with a microplate reader (ELX808; BioTek). The data were normalised to the control group.

Statistical analysis

All statistical analyses were performed using GraphPad Prism v.7.0 (GraphPad Software, La Jolla, CA, USA). All experiments were done at least three times. The data is shown as the mean ± SD. When comparing two groups, statistical significance was established using a two-tailed Student’s t-test, and when comparing more than two groups, one-way ANOVA followed by Tukey’s post-hoc test was used. A p-value <0.05 was used to signify statistical significance.

Similarity-based metric identifies osteogenic differentiation-related genes

To identify osteogenic differentiation-related genes, we performed similarity-based analysis to measure the topological relevance of genes involved in the human PPI network. It can be divided mainly into three steps: data collection, similarity measurement, and osteo-specific scoring. In the first step, human PPI network involving 19,344 genes is included in our analysis which then simplified by removing multiple edges and self-loops (Fig. 1B). Additionally, osteogenic differentiation-related genes are retrieved from published databases including GO, KEGG, BIOCARTA, PID, REACTOME, and WikiPathways (Fig. 1C). These genes are then used to extract osteogenesis-related features and grade genes into different levels. In the second step, we proposed two indexes Jaccard similarity and Sorensen-Dice similarity analysis to evaluate the topological similarity of genes in human PPI network (Adamic & Adar, 2003) (Fig. 1A). These two similarity indexes mainly concentrate on the overlap of neighbors of vertex A and B and measure the proportion of co-neighbors to evaluate the affinity of genes. After evaluating the similarity of genes using Jaccard and Sorensen-Dice similarity, they are then integrated into a final similarity matrix that presents the topological relevance of genes in human PPI network. In the third step, we mainly focus on these osteogenic differentiation-related genes and check whether some genes are significantly associated with these osteo-specific genes. We calculated the sum of similarity scores with these osteogenic differentiation-related genes for each gene which is then used to rank genes into different orders (Fig. 1D). This ranked gene list with osteo-specific score has potential in identification of osteogenic differentiation-related genes. We note that these genes with larger osteo-specific scores are more likely to participate in the regulation of osteogenic differentiation-related activity than other genes.

Reliable of these candidate genes in regulating osteogenic differentiation of hBMSCs

A total of 337 osteogenic differentiation-related genes are involved in our analysis to grade genes into different levels. All of these genes are retrieved from osteogenic differentiation-related pathways. But whether these genes contribute equally to the osteo prioritization of genes is unclear. To identify these significant determinants of osteogenic differentiation-related genes, we selected three pathways involving osteoblast-related transcription factors, osteoblast differentiation, and RUNX2 regulation of osteoblast differentiation for investigation. We found that RUNX2, CTNNB1, JUN, and FOS are significant regulators of transcription factors (Fig. 2A). Analysis of the osteoblast differentiation pathway shows that BGLAP, SPP1, HOX, HGF, IGF1, BMP2, BMP4, SMAD3, and NOTCH1 are significant regulators (Fig. 2B). Analysis of RUNX2 regulation of osteoblast differentiation pathway shows that RUNX2, BGLAP, COL1A1, SRC, MAPK3, AR, MAPK1, and ABL1 are significant regulators (Fig. 2C). In summary, these analyses suggest that these genes are significant regulators of osteogenic differentiation and play significant roles in the prioritization and determination of osteogenic differentiation-related genes.

To check whether these genes identified by similarity metrics can capture osteogenic features and are associated with osteogenic differentiation of hBMSCs, we investigated them in published literature. We found that almost all of these osteogenic differentiation-related genes lie in the upper gene list, suggesting that this ranked gene list can correctly prioritize and identify ossification-related genes. After removing these osteogenic differentiation-related genes, SMAD2, KLF4, DKK1, FGF13, NES, CCND1, SMAD7, IGF1R, TGFBR1, and KDR rank as the top 10 genes. Literature exploration shows that almost all of these genes have been reported to be associated with osteogenic differentiation. Activation of SMAD2 signaling pathway has been found to be associated with osteogenic differentiation (Zheng et al., 2020). KLF4, which encodes a protein that belongs to the Kruppel family of transcription factors, has been discovered to be a novel transcription factor of osteoblast differentiation (Yu et al., 2021). It has been found to regulate diverse cellular processes including cell proliferation, differentiation, and growth (Ghaleb & Yang, 2017). CCND1, which belongs to the highly conserved cyclin family, has been widely studied in various literature. Its downregulation has been observed to repress osteogenic differentiation and proliferation (Wang & Cai, 2020). Targeting Smad7 has been found to repress osteogenic differentiation of BMSCs (Fang et al., 2019a). The regulatory axis of GR/let-7f-5p/TGFBR1 has been discovered to be important for Dex-inhibited osteoblast differentiation (Shen et al., 2019). The IGF1R/PI3K/Akt signaling pathway has been shown to play a pivotal role in regulating osteogenesis (Fang et al., 2019b). In conclusion, these evidence suggests that our similarity-based algorithm is reliable in the prioritization of osteogenic differentiation-related genes.

Functional enrichment analysis of these candidate genes

To check whether these candidate genes identified by similarity-based metrics are enriched in some specific biological functions and states, we performed functional enrichment analysis of these top-ranked genes. The top 50 genes are selected for investigation of functional enrichment analysis. We performed functional enrichment analysis of these genes in GO and KEGG databases separately. GO enrichment analysis of these candidate genes shows that they are mainly enriched in osteogenic differentiation-related function including ossification and osteoblast differentiation (Fig. 3A). Additionally, other pathways such as mesenchyme development and mesenchymal cell differentiation are also presented by these candidate genes. Wnt signaling pathway, which has been widely studied and discovered to be associated with osteogenic differentiation, is observed to be enriched by these candidate genes. KEGG enrichment analysis of these genes shows that they are mainly enriched in cancer-related pathways including gastric cancer, hepatocellular carcinoma, breast cancer, and basal cell carcinoma (Fig. 3B). Apart from these pathways, the signaling pathway that regulates pluripotency of stem cells is also enriched in our analysis, suggesting significant value of these genes in regulating stemness. Other pathways such as hippo signaling pathway and wnt signaling pathway that have previously been reported to be associated osteogenic differentiation are also presented by these candidate genes. The top six enriched GO terms of these genes, which include cell fate commitment, epithelial to mesenchymal transition, mesenchymal cell differentiation, mesenchyme development, osteoblast differentiation, and regulation of animal organ morphogenesis, are presented in Fig. 3C. These results collectively demonstrate that these candidate genes can indeed capture osteogenic differentiation-related features of hBSMCs and confirm again that our similarity-based methodology is reliable in the identification of osteogenic genes.

FOXA1 is a significant osteogenic differentiation-related transcription factor

Among these top 100 candidate genes, we found that some of them are transcription factors such as BMI1, FOXA1, SRF, MYC, FOXC2, CREBBP, and ZEB2. Literature mining of these transcription factors shows that most of these genes have been reported to be associated with osteogenic differentiation. FOXA1, which encodes a member of the forkhead class of DNA-binding proteins, has no literature reporting its interaction with osteogenic differentiation. Due to the significant rank of FOXA1 in the gene list, it is necessary and significant to investigate the potential value of FOXA1 in osteogenic differentiation.

To check whether FOXA1 can participate in the regulation of osteogenic differentiation-related activity, we selected four well-known pathways for exploration which include the hedgehog signaling pathway, BMP signaling pathway, ERK pathway, and Wnt signaling pathway that have widely been reported to be related to osteogenic differentiation. Permutation test is first performed to check whether FOXA1 is significantly associated with these four pathways. We found that these four pathways have significantly higher similarity with FOXA1 than random noise (Fig. 4A), suggesting that FOXA1 is a significant regulator in mediating these osteogenic differentiation-related pathways.

We hypothesize that some genes involved in these pathways are more likely to be regulated by transcription factor FOXA1. In that way, FOXA1 can participate in the process of osteogenic differentiation-related activity by indirectly mediating these pathways. To identify target genes that are significantly regulated by transcription factor FOXA1, we analyzed the similarity score of these pathway genes with FOXA1. We found that GATA4, DKK1, SMAD2, and FOXD1 are significantly associated with FOXA1 in the BMP pathway with higher similarity scores than other genes (Fig. 4B), suggesting that these genes can be potentially regulated by FOXA1. Analysis of the Wnt signaling pathway shows that CCND2, DKK1, TCF7, HDAC2, NCOR2, and MYC are potential targets of FOXA1 (Fig. 4C). Analysis of the other two pathways hedgehog signaling pathway and erk pathway shows that VEGFA, ETS2, TLE1, TLE3, and NKX6-1 are potential targets in the hedgehog pathway, and STAT3, EGFR, IGF1R, PDGFRA, and MAP2P1 are potential targets in the ERK pathway (Figs. 4D–4E). These analyses suggest that FOXA1 participates in the regulation of osteogenic differentiation by indirectly regulating these target genes of ossification-related pathways.

FOXA1 knockdown enhanced alkaline phosphatase (ALP) activity and calcium deposit formation whereas FOXA1 overexpression decreased ALP activity and calcium deposit formation

ALP activity is a characteristic of early osteogenesis. On the 5 day of osteogenic differentiation, ALP activity was evaluated. ALP activity was significantly higher in the FOXA1-KD group compared to the KD-NC group (p < 0.05), and ALP staining results were similar (Figs. 5A and 5B). In comparison to the FOXA1 OE-NC group, the FOXA1 OE group had decreased ALP activity (P < 0.05) (Figs. 5A and 5B).

The calcium deposits were investigated using alizarin red staining (ARS). On day 16, the FOXA1 knockdown group had more calcium deposits than the KD-NC group, whereas the FOXA1 OE group had less calcium deposits than the OE-NC group. Quantification analysis yielded comparable results (Figs. 5A and 5B).