Identification of a potential diagnostic signature for postmenopausal osteoporosis via transcriptome analysis.

Purpose: We aimed to establish the transcriptome diagnostic signature of postmenopausal osteoporosis (PMOP) to identify diagnostic biomarkers and score patient risk to prevent and treat PMOP.
Methods: Peripheral blood mononuclear cell (PBMC) expression data from PMOP patients were retrieved from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) were screened using the "limma" package. The "WGCNA" package was used for a weighted gene co-expression network analysis to identify the gene modules associated with bone mineral density (BMD). Least absolute shrinkage and selection operator (LASSO) regression was used to construct a diagnostic signature, and its predictive ability was verified in the discovery cohort. The diagnostic values of potential biomarkers were evaluated by receiver operating characteristic curve (ROC) and coefficient analysis. Network pharmacology was used to predict the candidate therapeutic molecules. PBMCs from 14 postmenopausal women with normal BMD and 14 with low BMD were collected, and RNA was extracted for RT-qPCR validation.
Results: We screened 2420 differentially expressed genes (DEGs) from the pilot cohort, and WGCNA showed that the blue module was most closely related to BMD. Based on the genes in the blue module, we constructed a diagnostic signature with 15 genes, and its ability to predict the risk of osteoporosis was verified in the discovery cohort. RT-qPCR verified the expression of potential biomarkers and showed a strong correlation with BMD. The functional annotation results of the DEGs showed that the diagnostic signature might affect the occurrence and development of PMOP through multiple biological pathways. In addition, 5 candidate molecules related to diagnostic signatures were screened out.
Conclusion: Our diagnostic signature can effectively predict the risk of PMOP, with potential application for clinical decisions and drug candidate selection.

Introduction

Osteoporosis is the most common systemic bone disease in postmenopausal women (PMOP). This condition is characterized by decreased bone mineral density (BMD) and destruction of the bone tissue microstructure, resulting in decreased bone strength and increased fracture risk (Yao et al., 2019). The latest data show that the prevalence of osteoporosis in women over 50 years old in China is 29.1%, and the total number of patients is approximately 49 million. With the increase in the aging of society, by 2050, there will be 5.99 million cases of osteoporosis-related fractures in China every year (Cheng et al., 2021). Hip fracture is a common osteoporotic fracture in older postmenopausal women. Approximately 33% of women aged 90 years will experience a hip fracture, and these patients are older and have multiple comorbidities, leading to a poor prognosis (Ensrud et al., 2019). Osteoporosis has become a major public health problem worldwide, placing a heavy financial burden on patients and healthcare systems. Osteoporosis includes primary and secondary osteoporosis. The occurrence of primary osteoporosis is related to estrogen deficiency in women, decreased testosterone levels in men, and changes in hormone levels such as parathyroid hormone and calcitonin. Secondary osteoporosis is mainly caused by endocrine and metabolic diseases, connective tissue diseases, kidney diseases, digestive tract diseases, and drugs.

Most patients have no apparent symptoms in the early stage of osteoporosis, and the main clinical manifestations include pain, reduced height, limited activities, stooped posture, and respiratory system involvement. Most people lack awareness of osteoporosis and fail to detect the symptoms. Even when a brittle fracture occurs, there is no clear history of trauma or only a slight history of trauma (Glaser and Kaplan, 1976). Osteoporosis can be characterized by sparse trabecular bone and decreased BMD on imaging, but these imaging manifestations are affected by subjective factors and are not sensitive to early bone loss. Dual-energy X-ray absorptiometry (DXA) testing of BMD is strongly recommended by the World Health Organization (WHO). However, DXA cannot detect early bone loss, and studies of many clinical cases have suggested that DXA does not accurately assess the severity of osteoporosis and the risk of fracture. For example, in one report, the individuals who had the highest risk of future fractures among those who had bone density tests were rarely diagnosed with osteoporosis (T < -2.5) but were often diagnosed with reduced bone mass (−2.5 < T < -1) (Siris et al., 2014).

Therefore, transcriptome analysis may be helpful for the early diagnosis and prevention of osteoporosis. Based on the Gene Expression Omnibus (GEO) database, we used WGCNA to identify the top gene modules related to osteoporosis and conducted least absolute shrinkage and selection operator (LASSO) regression to build a diagnostic signature composed of 15 genes. The two most representative genes, METTL4 and RAB2A, were subjected to RT–qPCR to verify the correlation between their expression and BMD. In addition, we explored the potential molecular mechanism of osteoporosis, which will contribute to the early diagnosis, treatment, and prevention of this disease. To discover novel osteoporosis drugs from our research, we explored the molecular targets of the diagnostic signature through network pharmacology.

Materials and methods

Microarray data collection and processing

The datasets were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013), and the pilot cohort and discovery cohort were analyzed on the Affymetrix Human Genome U133A Array platform. As the pilot cohort, GSE56815 included 80 Caucasian females, 40 of whom had high BMD and 40 of whom had low BMD. In addition, we selected 20 postmenopausal subjects from the GSE13850 dataset as a discovery cohort: 10 women with high BMD and 10 women with low BMD. Raw data were read through the “affy” package (Gautier et al., 2004), and the RMA algorithm was used for background correction and data normalization. To verify the purity of peripheral blood mononuclear cells (PBMCs) in the pilot cohort, we used the “Cibersort” package.

Identification of differentially expressed genes (DEGs)

Patients are classified into high- and low-BMD groups based on clinical information provided by uploaders in the pilot cohort. The “limma” package (Ritchie et al., 2015) was used to compare samples from the high- and low-BMD groups. Genes with adjusted p-value < 0.05 were defined as DEGs, DEGs with log2FC > 0 were defined as up-regulated DEGs, while those with log2FC < 0 were down-regulated DEGs. The “ggplot2” package (Ito and Murphy, 2013) was used to visualize the DEGs.

Weighted gene coexpression network analysis

WGCNA is an algorithm to mine module information from high-throughput data. In this method, the module is defined as genes with similar expression trends. If these genes always have similar expression changes in a physiological or pathological process, it is reasonable to believe that they are functionally related and can be defined as a module. In this study, the “WGCNA” package (Langfelder and Horvath, 2008) was used to construct the weighted adjacency matrix by selecting appropriate thresholds, and the weighted adjacency matrix was transformed into a topological overlap matrix (TOM). The hierarchical clustering method was used to cluster the TOM matrix, and the dynamic tree-cutting algorithm was adopted to divide modules, each corresponding to a color, merge similar modules to find the module with the highest correlation with BMD, and extract the most significant genes associated with BMD in the module.

Construction and validation of a genetic diagnostic signature of osteoporosis

LASSO logistic regression was used to reduce the dimensionality of genes in the BMD association module to construct a genetic diagnostic signature, which was generated by using the “glmnet” package (Friedman et al., 2010). Then, the signature was used to calculate the risk score of each patient. The corresponding coefficient of the gene weighted the expression values of these genes in each patient, and then, the weighted expression values were added to obtain the risk score of the patient, which was calculated as follows: where n is the number of genes included in the signature, is the expression value of this gene of the patient, and is the coefficient of this gene in the signature. Finally, patients were classified into a high-risk group or a low-risk group according to the median value of the risk score.

Next, the “pROC” package (Robin et al., 2011) was used to draw the receiver operating characteristic curve (ROC) and determine the area under the curve (AUC). If the AUC was >0.8, the diagnostic effect was defined as good. In addition, we obtained another independent dataset, GSE13850, as the discovery cohort, applied the genetic diagnostic signature to the discovery cohort and defined its diagnostic effect in the discovery cohort according to the AUC.

Functional annotation analysis

For determination of the potential biological pathways of osteoporosis, “clusterProfiler” package (Yu et al., 2012) was used to perform Kyoto Encyclopedia of Genes (KEGG), Gene Ontology (GO), Disease Ontology (DO) on the up- and down-regulated DEGs, the results for which adjust p-value<0.05 were considered statistically significant. Gene Set Enrichment Analysis (GSEA) on the up- and down-regulated DEGs was performed by “fgsea” package.

Gene-miRNA interaction analysis, molecular docking, and network visualization

TargetScan (McGeary et al., 2019) was used to predict the interacting miRNAs of the genes in the genetic diagnostic signature and explore whether potential miRNAs are involved in the process of osteoporosis by targeting genes. Subsequently, miRPath was used for pathway enrichment analysis of miRNAs (Vlachos et al., 2012).

Firstly, we screened potential small molecule substances and drugs from the CTD database (Davis et al., 2021) and DGIdb database. Cytoscape was used to construct the miRNA-gene-molecule network. Secondly, we downloaded the protein 3D structures from the PDB database (Burley et al., 2021) and used PyMOL software to remove the water molecules of the protein. Then, we collected the 3D structures of these potential small-molecule substances and drugs from PubChem. Next, the SwissDock (Grosdidier et al., 2011) database was used for protein-molecule docking. Finally, PyMOL software was used to modify the docking results.

RNA extraction and quantitative real-time polymerase chain reaction

Fourteen women with PMOP and 14 healthy postmenopausal-matched women controls were selected, and PBMCs were obtained from them. The characteristics of the patients are summarized in Table 1. This study was approved by the Ethics Committee of the First Affiliated Hospital of Jinan University. Both patients and controls provided written informed consent.

TABLE 1

Patient characteristics (n = 28).

	Mean ± SD (Range)
Characteristic	High bone density	Low bone density
Patients (n)	14	14
Age (years)	61.9 ± 9.5	66.2 ± 8.4
Height (cm)	156.8 ± 6,7	158.6 ± 5.6
Weight (kg)	58.4 ± 6.3	61.1 ± 7.6
T-score	0.7 ± 0.9	-2.4 ± 0.8
Menstrual condition
Menopause	14	14
Premenopausal	0	0
Smoking status
Smoking	0	0
No history of smoking	14	14
Surgery situation
Have vertebroplasty	0	0
No vertebroplasty	14	14

PBMCs were extracted with Histopaque-1077 (Sigma, United States). According to the manufacturer’s protocol, the total RNA of PBMCs from all samples was extracted using an EZ-Press RNA Purification Kit (EZbioscience, United States). cDNA was obtained by reverse transcription using a PrimeScript RT Kit (TaKaRa, Japan). Based on the SYBR Green method (ChamQ Universal SYBR qPCR Master Mix, Vazyme Biotech, China), The CFX96 Real-Time PCR System (Bio-Rad, United States) was used for RT-qPCR detection. The mRNA-specific primer sequences are shown in Table 2. After the expression level of GAPDH was used for normalization, the relative expression level of mRNA was determined.

TABLE 2

mRNA-specific primer sequences.

Gene	Primer sequence	Tm
METTL4	F: GCT​GTT​CAT​AAA​GAA​TGC​CAG​CAA	57
	R: CAG​CTC​CCT​GAT​CTT​TGT​ATG​GT	56
RAB2A	F: TCC​ATC​ACA​AGG​TCG​TAT​TAC​AGA	55
	R: TGG​TTG​AAT​GTA​TCT​CTC​CGT​GTA	55
GAPDH	F: ACAGTTGCCATGTAGACC	56
	R: TTTTTGGTTGAGCACAGG	60

Statistical analysis

Statistical analyses were conducted using RStudio version 1.3.1093 (RStudio Inc.) and GraphPad Prism 8 (GraphPad Software, Inc.). All data are expressed as the mean ± SD. A paired difference test between Low BMD samples and normal control samples in two groups by the “limma” package was used to determine the DEGs. Analysis with one-way ANOVA followed by the Student–Newman–Keuls multiple comparison test was used for the comparison of three or more experimental groups. For qPCR data, Student’s t test was used for analysis.

Results

Quality control of microarray data and DEG screening

First, cell purity analysis was performed on the pilot cohort (GSE56815), and the results are shown in Figure 1A, indicating that PBMCs accounted for the majority of the pilot cohort, consistent with the description in NCBI. DEG analysis was performed on the pilot cohort, and patients were divided into two groups according to BMD. DEGs were identified by “limma” package analysis, and 2,420 DEGs were screened out, as shown in Figure 1B.

FIGURE 1

(A) Cell abundance of the pilot cohort. (B) Volcano plot of DEGs between the individuals with a high BMD and a low BMD in the pilot cohort.

Weighted coexpression network and identification of bone mineral density-related modules

In this study, the “WGCNA” package was used to construct a weighted gene coexpression network in the pilot cohort. According to several iterations, β = 6 was selected as the optimal soft threshold to construct a scale-free network (Figure 2A). After exclusion of the MEgrey module, which contains all genes not involved in clustering, a total of 7 modules were identified. The interaction between modules was analyzed, and the heatmap showed that the gene expression of each module was relatively independent (Figure 2B). Then, correlation analysis between these modules and BMD was carried out, and the results showed that the MEblue module, which consisted of 396 genes, had the highest correlation with BMD (cor = 0.51, p < 0.001) (Figure 2C).

FIGURE 2

(A) Fitting index of the scale-free topology module under different soft thresholds (left) and network connectivity under different soft thresholds (right). (B) Cluster diagram of genes and the corresponding gene modules. (C) Correlation between module eigengenes and BMD. (D) Selection of the optimal parameter (λ) of LASSO regression through cross-validation. (E) LASSO coefficient profiles of the 15 genes that comprise the diagnostic signature selected by λ.

Construction and testing of the diagnostic signature

The LASSO algorithm was used to determine λ = 0.09 (Figure 2D), and a diagnostic signature consisting of 15 genes (Figure 2E) was established. The specific gene composition and coefficient of each gene are shown in Table 3. The diagnostic signature calculated the patients’ risk scores and divided them into high-risk and low-risk groups (Figure 3A). Principal component analysis (PCA) showed that risk scores could categorize patients with different BMDs in the pilot cohort into two groups (Figure 3B). Moreover, verification was carried out in this study, and the results showed that the AUC value of the diagnostic signature was 0.993 in the pilot cohort and 0.980 in the discovery cohort (Figure 3C). The expression patterns of the 15 genes that constituted the signature (Figure 3D) and all the DEGs (Figure 3E) in the high-risk and low-risk groups are shown. Heatmap results showed that the gene expression patterns of the high-risk and low-risk groups were different, especially those of the 15 genes that constituted the signature. All the above results showed that this signature has an excellent ability to predict the risk of osteoporosis. In addition, we explored the interactions of these 15 genes (Figure 3F).

TABLE 3

Genes and their coefficients that constitute the diagnostic signature.

Gene	Coefficient
RAB2A	−0.69559
VSIG4	0.005073
ADAM7	0.014223
AMBP	0.018658
PAFAH1B2	0.031496
AOC3	0.050969
KLK3	0.088715
KRTAP1.3	0.100,972
LPO	0.13443
SLC41A3	0.144,549
NKX3.1	0.208,241
GLT8D2	0.231,434
LAMB1	0.338,636
SEC14L1P1	0.394,242
METTL4	1.089913

FIGURE 3

(A) Risk plot and (B) PCA of the pilot cohort. (C) The heatmap shows the different expression patterns of 15 genes. (D) ROC curves of the pilot cohort and the validation cohort. (E) The heatmap shows the different expression patterns between the high- and low-risk groups. (F) Correlations among the 15 genes.

Functional annotation and Gene Set Enrichment Analysis

KEGG functional annotation analysis (Figure 4A) showed that cytokine receptor interaction, neural activity ligand–receptor interaction, Rap1 signaling pathway, autoimmune thyroid disease, natural killer cell-mediated cytotoxicity, PI3K-Akt signaling pathway, gap junction, calcium signaling pathway, and oxytocin signaling pathway were up-regulated. Combining the above results with previous studies, we found that the calcium activator calcimycin can activate the RAF/MEK/ERK pathway through Ras (Li et al., 2005), and increased calcium concentrations have also been shown to modulate Ras-dependent Raf1 activation (Yoshiki et al., 2010). Moreover, lactoferrin-induced PI3K-Akt pathway activation and Ras phosphorylation can promote osteoblast proliferation (Hou et al., 2015). DO analysis (Figure 4B) showed that gynecological and aging diseases were up-regulated, such as female reproductive system disease, ovarian disease, bone remodeling disease, osteoporosis, and bone resorption diseases. GO analysis (Figure 4C) showed that ion channel complex activity, bone development, bone morphogenesis were up-regulated, and odontogenesis, GTPase activity, GDP binding were down-regulated.

FIGURE 4

(A) Result of KEGG enrichment analysis. (B) Result of DO analysis. (C) Result of GO analysis. (D) Results of the GSEA.

According to the GSEA results (Figure 4D), we suggest that androgen response, late estrogen response, the P53 pathway, and TNF-α signaling via NF-κB are were down-regulated in the DEGs of osteoporosis.

miRNA interaction identification and candidate molecule prediction

Thirty-eight miRNAs were expected to interacted with the genes constituting the diagnostic signature and were used to construct the miRNA-Gene-Molecule network (Figure 5A). The result of KEGG (Figure 5B) enrichment analysis showed that miRNA pathway enrichment overlapped with DEG pathways, such as the Rap1 signaling pathway, suggesting that overlapping pathways play a potentially important role in the occurrence and development of PMOP. We selected the molecules connected to at least two genes as the candidate molecules. Five candidate molecules were screened out, including bisphenol A (BPA), fulvestrant, bicalutamide, mifepristone, and valproic acid (VPA). RAF, an essential protein in the Rap1 signaling pathway, was selected for docking with the candidate molecules to explore their possible binding locations (Supplementary Figure S1).

FIGURE 5

(A) miRNA-gene-molecule interaction network. (B) KEGG enrichment analysis of the potential interacting miRNAs.

Validation by real-time polymerase chain reaction

To verify the authenticity of the diagnostic signature, we collected PBMCs from 14 postmenopausal healthy controls (Figure 6A) and 14 postmenopausal women with low BMD (Figure 6B) in this study. RNA was extracted for RT–qPCR to verify the diagnostic signature. The gene with the most significant coefficient had the strongest contribution to the risk score, and the gene with the largest AUC showed the strongest relationship to BMD. Interestingly, METTL4 and RAB2A had the largest positive and negative coefficients, respectively, and they also had the most significant AUC values (Table 4). Therefore, these two genes were further investigated, and the results were consistent with the bioinformatics results. The expression of METTL4 in the low BMD group was significantly higher than that in the normal BMD group (Figure 6C), while the expression of RAB2A in the normal BMD group was significantly higher than that in the low BMD group (Figure 6D), suggesting that these two genes play an essential role in the occurrence and development of osteoporosis.

FIGURE 6

DXA images of the lumbar spine L1-L4 of women with a normal BMD (A) and low BMD (B). The overall expression of METTL4 (C)and RAB2A (D) in PBMCs from 14 low BMD patients and healthy controls.

TABLE 4

Genes and their AUCs that constitute the diagnostic signature.

Gene	AUC
RAB2A	0.8638
VSIG4	0.7431
ADAM7	0.7219
AMBP	0.72
PAFAH1B2	0.725
AOC3	0.6444
KLK3	0.6444
KRTAP1.3	0.7375
LPO	0.7141
SLC41A3	0.7338
NKX3.1	0.7075
GLT8D2	0.7819
LAMB1	0.6481
SEC14L1P1	0.7038
METTL4	0.8306

Discussion

Osteoporosis is a systemic bone disease that mainly involves decreased bone mass and increased bone brittleness caused by degeneration of the bone tissue microstructure, resulting in susceptibility to fracture (Wang et al., 2020). DXA is widely used in clinical practice as a diagnostic standard for osteoporosis. However, only a small number of postmenopausal women are tested for BMD, and many of them have already suffered from brittle fractures when BMD is found to be reduced (Siris et al., 2014). BMD examination results may be affected by body weight, lumbar curvature, osteophytes, vertebral fractures, and vascular calcification. Studies have shown that vertebral osteoarthropathy and aortic calcification can cause false BMD increases and decrease diagnostic sensitivity (Orwoll et al., 1990; Frohn et al., 1991).

Relying solely on BMD testing will lead to the failure of clinicians and patients to correctly evaluate the severity of osteoporosis, thus affecting clinical treatment strategies. In addition, the radiation of DXA is very low, with each scan receiving only 1/60 to 1/10 of the radiation dose of a conventional X-ray. However, repeated exposure to ionizing radiation over a long period can have long-term health effects, including cancer (Hill and Einstein, 2016; Howard et al., 2020), suggesting that frequent imaging detection is still not advisable.

Due to the difficulty in the extraction and separation of osteoblasts and osteoclasts, on the contrary, the isolation and extraction technology of PBMCs has become increasingly mature, and its extraction and purification rate can reach more than 90%. PBMCs are the most likely precursors of osteoclasts, especially in adult peripheral bone, and are the only precursors of osteoclasts (Fujikawa et al., 1996). Secondly, PBMCs can secrete cytokines such as IL-1B, IL-6, and TNF-α, which play an essential role in osteoclast differentiation, activation, and apoptosis (Custer and Ahlfeldt, 1932). The decrease of PBMCs cytokines is the primary mechanism by which sex hormones inhibit osteoclast formation and bone resorption. Therefore, PBMCs are widely used as an ideal cell model for osteoporosis study (Zhou et al., 2015).

Previous studies have shown that osteoporosis is a polygenic disease, and genetic factors play an essential role in the occurrence and development of this condition. However, there are few relevant studies at present, and few candidate genes have been identified (Andrew and Macgregor, 2004; Cai et al., 2018). Thus, we used public databases to construct a transcriptome-based diagnostic signature, which is expected to be applied in clinical practice in the future to assess the risk of osteoporosis in the potential population and help doctors diagnose this disease. In this study, the GEO dataset GSE56815 was first used to identify DEGs between patients with osteoporosis and healthy controls, and then WGCNA was performed on these DEGs. The results showed that the blue module was most closely related to osteoporosis. Then, the genes in the blue module were extracted, and a diagnostic signature composed of 15 genes was constructed through LASSO regression. Subsequently, the independent dataset GSE13850 was used to verify the classification ability of the signature. We found that the AUC values of the signature in the pilot cohort and the discovery cohort were 0.993 and 0.920, respectively, showing good prediction of the risk of osteoporosis in different patients.

To further explore the potential pathogenesis of osteoporosis, we conducted the functional enrichment analysis of DEGs between the high- and low-BMD groups and interacting miRNAs of the diagnostic signature. We found that both were involved with the Rap1 signaling pathway, suggesting that, compared with other mechanisms, the Rap1 signaling pathway may be more closely related to our diagnostic signature. The calcium activator calcimycin can activate the RAF-MEK-ERK pathway through the RAS signaling pathway (Li et al., 2005). An increased calcium concentration has also been shown to modulate RAS-dependent RAF1 activation (Yoshiki et al., 2010), and lactoferrin-induced PI3K-Akt pathway activation and Ras phosphorylation can promote osteoblast proliferation (Hou et al., 2015). According to previous studies, osteoclasts are specialized macrophage/monocyte lineage-derived cells that resorb bone, and neurofibromatosis type I (NF1) haploinsufficient osteoclasts have abnormal Ras-dependent bone resorption. (Yan et al., 2008).

Our study found that the two genes with the most significant positive and negative coefficients in the diagnostic signature, and , were also the two genes with the most significant AUC values among the 15 genes. Therefore, we selected these two genes as representative molecules that were validated by RT–qPCR. We collected PBMCs from 14 patients with normal BMD and 14 patients with low BMD for RNA extraction. The results showed that expression was significantly lower in patients with normal BMD than in patients with low BMD, while expression was significantly higher in healthy controls than in patients with low BMD, which was consistent with our bioinformatics results and the results of other studies. is a member of the Ras gene family. When it binds to GTP, the Ras protein can phosphorylate and activate downstream proteins, thus regulating the proliferation and differentiation of osteoblasts (Ge et al., 2007). The Ras gene family regulates anterograde transport from the endoplasmic reticulum to the Golgi complex, inhibits the proliferation and differentiation of osteoprogenitor MC3T3-E1 cells and promotes apoptosis by reducing the membrane transport process. During this process, the protein expression of was inhibited (Hong et al., 2011). METTL4, as an m6A methyltransferase, can lead to increased m6A modification. As a new extranuclear marker, m6A modification is involved in bone development and metabolism and plays an essential role in osteoporosis. A previous study showed that both m6A levels and methyltransferase expression are increased during osteoclast differentiation, and methyltransferase knockdown led to increased osteoclast volume and decreased bone resorption capacity (Li et al., 2020). Another study also showed that inhibition of mRNA methyltransferase reversed osteoclast differentiation and bone resorption (Wang et al., 2021).

Another vital function of diagnostic signatures is to provide evidence for candidate drugs. RAF is a crucial component of the Rap1 signaling pathway, and RAF has also been shown to affect the proliferation and function of osteoblasts (Meng et al., 2015). Therefore, it is necessary to develop effective new osteoporosis drugs that target RAF. This study screened out five medications with high affinity for the diagnostic signature: BPA, fulvestrant, mifepristone, bicalutamide, and VPA. Although the specific mechanisms of action of these small compounds remain to be further elucidated, our results suggest that they have therapeutic potential for osteoporosis, especially in patients with PMOP.

We speculated that using transcriptome analysis to detect gene expression might be a complementary method for the diagnosis of osteoporosis, and the genes can be used as biomarkers to evaluate the effect of osteoporosis treatment to avoid frequent radiation examinations in patients.