Introduction
Because of demographic shifts such as global population aging and increased life expectancy, osteoporosis has become a major public health concern worldwide. The incidence is significantly higher in women than in men, and the number of affected individuals is projected to double between 2010 and 2040 [ 1 ]. Osteoporosis is a metabolic bone disease characterised by reduced bone strength and increased fracture risk [ 2 ], often leading to a significant decline in quality of life. It is estimated that approximately 20% of individuals aged 50 and older who sustain a hip fracture will die within 1 year, and nearly 50% will experience long‐term disability [ 3 ]. Therefore, improving early diagnosis and developing more effective therapeutic strategies are critical.
The pathogenesis of osteoporosis is multifactorial, involving age, sex, genetic susceptibility, medication use and comorbid conditions such as rheumatoid arthritis [ 4 ]. Although current treatments–such as bisphosphonates and monoclonal antibodies–are clinically effective, they are associated with limitations and adverse effects, including osteonecrosis of the jaw and atypical femoral fractures [ 5 ]. Thus, identifying novel molecular contributors and therapeutic targets is essential to advancing treatment options and minimising complications.
Mendelian randomisation (MR) is a genetic epidemiological approach that uses genetic variants as instrumental variables to assess potential causal relationships between exposures and disease outcomes [ 6 ]. Drug target MR, an extension of this method, evaluates variants in or near druggable genes, thereby mimicking the biological effects of pharmacological interventions [ 7 ]. This approach has been applied to assess the efficacy and safety of medications such as angiotensin‐converting enzyme inhibitors [ 8 , 9 , 10 ], and is increasingly used in drug discovery and repurposing efforts [ 9 , 11 ].
Despite its value, conventional MR analyses are largely limited to gene–phenotype associations and often fail to provide mechanistic insight into downstream transcriptional regulation or cell‐type specificity. Furthermore, integrative frameworks combining MR with transcriptomic and single‐cell analyses remain rare. To address these gaps, we developed a comprehensive ‘MR–Transcriptome–Single Cell’ approach that (1) identifies candidate drug targets using FinnGen and EBI GWAS data, (2) maps biological pathways—such as bile acid signalling—regulated by FMO4 through gene set enrichment analysis (GSEA) and (3) locates FMO4 expression specifically in osteoblast precursor cells using single‐cell RNA sequencing data. This integrative strategy aims to enhance the identification of osteoporosis‐related targets and support future therapeutic research.
Results
Mendelian Randomisation Analysis of the Training Set
We obtained druggable genes from the literature [ 12 ]. We obtained osteoporosis‐related GWAS data from the FinnGen biobank to pinpoint key druggable genes affecting osteoporosis. The summary statistics for 399,054 samples (Controls: 391,037; Cases: 8017) were used, with the outcome ID for osteoporosis being finngen_R10_M13_OSTEOPOROSIS. Data for exposure and outcome were extracted using the extract_instruments and extract_outcome_data functions. MR analysis revealed 85 genes linked to eQTL‐positive outcomes (Figure 1 , IVW p < 0.05). The genes TUBB4B, PLA2G2C, C3AR1, PREP, NOG, H3F3B, ALDH8A1, ITGA3, CATSPER1, COL17A1, PRTN3, KALRN, SCGB1C1, TLR10, TNFRSF13B, FMO4, CTLA4, KL, GGH, CCL3L1, SPON2, CRTAM, TPBG, NOV, XYLT1, PNP, CD302, ALDH5A1, IL15, SRPK1, ITPR3, CES1, KCNK6, MARK3, PINK1, SLC7A8, ATF1, CORIN, TNFSF4, SLC18A1, APAF1, POGLUT1, OPN3, ITGAX, DCPS, GYPE, and TRPM6 are linked to a reduced risk of osteoporosis, whereas PVRL2, OAS1, BTNL8, GSTM1, SERPINB1, TUBB2A, CCR1, RTN4, IL12RB2, CCL23, STK11, SCARB1, CTNNB1, MAP2K5, EPOR, PARP1, CTSF, SDC2, SLAMF1, VEGFA, MAP3K1, MYLK4, MAP4K2, HMGCR, ICOS, APEX1, CCR9, PSMA4, TUBB, FCGR1A, CYP21A2, MGAT4A, AQP1, PTAFR, UPP1, SUCNR1, KARS and ANGPT2 are associated with an increased risk of osteoporosis. Additional sensitivity analysis of the 85 causal relationships indicated that excluding any single SNP did not notably impact the overall error line. This indicates that the 85 causal relationships selected in this study are robust.
MR Analysis and Colocalisation Analysis of the Validation Set
The validation set osteoporosis‐related GWAS summary data were obtained from 667,227 samples (controls: 648,913; cases: 18,314), with the outcome ID being GCST90018887. The data for exposure and outcome were sequentially extracted using the extract_instruments and extract_outcome_data functions.
MR analysis revealed six genes linked to eQTL‐positive outcomes (Figure 2A–F , IVW p ‐value < 0.05). FMO4 and NOV were associated with a decreased risk of osteoporosis, whereas MAP2K5, PSMA4, SCARB1 and VEGFA were linked to an increased risk. Additional sensitivity analysis of the six causal relationships demonstrated their robustness, as excluding any individual SNP did not notably alter the overall error line (Figure 3A–F ). The six genes were analyzed for eQTL‐GWAS co‐location. FMO4, PSMA4 and VEGFA showed high probabilities (> 0.75; Figure 4A–C ) and were selected as key genes for further analysis.
GSEA
We conducted GSEA to investigate the signalling pathways associated with the three key genes. The results indicated that FMO4 was enriched in the mRNA surveillance, calcium signalling, and cell adhesion molecules pathways (Figure 5A,B ). PSMA4 was enriched in fatty acid metabolism, Rap1 signalling pathway and PI3K‐Akt signalling pathway (Figure 5C,D ). VEGFA was enriched in the IL‐17 signalling pathway, MAPK signalling pathway and lysosome pathway (Figure 5E,F ).
Gene Set Variation Analysis (GSVA)
GSVA indicated that the FMO4 high‐expression group was enriched in bile acid metabolism, hedgehog signalling and Wnt‐beta‐catenin signalling pathways (Figure 6A ). PSMA4 was enriched in DNA repair, TGF‐beta signalling and heme metabolism pathways (Figure 6B ). VEGFA was involved in the signalling pathways of Wnt‐beta‐catenin, IL‐6/JAK/STAT3 and IL‐2/STAT5 (Figure 6C ). The results indicate that pivotal genes might affect disease progression via these signalling pathways. Notably, the observed enrichment of the ‘bile acid metabolism’ pathway may reflect broader systemic metabolic effects rather than bone‐specific alterations. This caveat should be considered when interpreting its relevance to osteoporosis pathophysiology.
Immune Infiltration Analysis
The local microenvironment, consisting of immune cells, extracellular matrix components, diverse growth factors, inflammatory mediators and specific physicochemical characteristics, exerts a profound influence on disease diagnosis and sensitivity to clinical treatment. In this study, we systematically analysed the relationship between key genes and immune cell infiltration within an osteoporosis dataset, aiming to further elucidate the potential molecular mechanisms through which these genes contribute to the progression of osteoporosis. Additionally, we characterised the proportional distribution of immune cell populations in each patient and investigated the correlations among different immune cell types (Figure 7A,B ). Moreover, statistically significant differences were identified between the two groups in terms of APC co‐stimulation, immature dendritic cells (iDCs), T follicular helper (Tfh) cells and T helper 1 (Th1) cells (Figure 7C ).
The Association Between Essential Genes and Immune Cells
We additionally investigated the association between essential genes and immune cells. FMO4 was found to be significantly positively correlated with CCR, CD8 + T cells, APC_co_inhibition and T helper cells (Figure 7D ). PSMA4 exhibited a notable positive correlation with Treg, MHC class I, and Type I IFN response, while demonstrating a significant negative correlation with T helper cells, CCR and APC_co_inhibition (Figure 7E ). VEGFA showed a significant positive correlation with macrophages and T helper cells, while exhibiting a negative correlation with cytolytic activity, DCs and inflammation‐promoting pathways (Figure 7F ).
Association Between Essential Genes and Immune Modulators
We examined the correlation between key genes and diverse immune factors, such as immunosuppressive and immunostimulatory factors, major histocompatibility complex (MHC), chemokines and receptors. These analyses indicate that key genes are significantly linked to immune cell infiltration levels and are essential in the immune microenvironment (Figure 8A–E ).
Transcriptional Regulation and miRNA Network Related to Key Genes
This study identified that the three key genes are regulated by multiple transcription factors and common mechanisms. Enrichment analysis was performed on these transcription factors using the cumulative recovery curve. Analysis of key genes for motif‐TF annotation and selection identified motif cisbp__M5572 as having the highest normalised enrichment score (NES = 5.76). Figure 9A,B illustrates the enriched motifs and associated transcription factors for the key genes identified in this study. Using the miRcode database, we predicted potential miRNAs for the three key genes, identifying 78 miRNAs and 123 mRNA‐miRNA pairs, which were visualised in Cytoscape (Figure 9C ).
Correlation Between Key Genes and Metabolic Pathways and Disease Progression Genes
Next, we analysed the correlation between the key genes and metabolic pathways. In the metabolic pathway heatmap, it was observed that the high‐expression group of FMO4 had a significantly higher score in the Drug metabolism relevant signatures category compared to the low‐expression group (Figure 10A ). We obtained osteoporosis‐related genes from the GeneCards database ( https://www.genecards.org/ ) and extracted the top 20 genes. A correlation analysis between the expression levels of key genes and disease‐related genes revealed significant associations. FMO4 showed a significant positive correlation with CYP19A1 (cor =0.37), whereas PSMA4 exhibited a significant negative correlation with ESR1 (cor = −0.526) (Figure 10B ).
Quality Control and Clustering of Single‐Cell Data
Following quality control, 7471 high‐quality cells were retained (Figure 11A,B ). These cells exhibited a median gene count > 200 and mitochondrial gene content within physiological ranges. DoubletFinder effectively removed putative doublets, ensuring dataset integrity. HVG analysis identified 2000 genes demonstrating significant expression variability across the cellular landscape (Figure 11C ). Dimensionality reduction through PCA (Figure 11D ) and batch effect correction using Harmony (Figure 11E ) facilitated downstream clustering.
Cell Type Annotation and Key Gene Expression Patterns
UMAP projection resolved 12 distinct clusters, which were annotated into seven major cell types: B cells, T cells, macrophages, neutrophils, vascular smooth muscle cells (VSMCs), plasmacytoid DCs and erythroid cells (Figure 12A,B ). Marker gene expression validated annotations, with canonical markers (e.g., CD79A for B cells, CD3D for T cells) displaying cell‐type‐specific enrichment (Figure 12C ). Key genes (FMO4, PSMA4 and VEGFA) exhibited differential expression patterns across cell types (Figures 12D,E ), reflecting lineage‐specific functional roles.
Drug Prediction Using the DGIdb Database
This study utilised the DGIdb database to examine potential drug interactions with key genes. Through DGIdb, 52 drugs were identified to interact with VEGFA, and 8 drugs were found to interact with PSMA4. These interactions may provide potential new targets for therapeutic development. Cytoscape was utilised to visualise the drug–gene interactions (Figure 13 ).
External Validation of the Identified Key Genes
To further validate the expression trends of the identified hub genes, we conducted cross‐dataset verification using the GSE35959 transcriptomic dataset, which comprises peripheral blood samples from osteoporosis patients and healthy controls. The analysis revealed that FMO4 and PSMA4 were consistently downregulated in osteoporosis patients, whereas VEGFA exhibited significantly increased expression levels (Figure 14 ). These findings align with the results from our primary analysis and provide additional support for the potential roles of these genes in the development and progression of osteoporosis.
Discussion
Osteoporosis is a systemic metabolic bone disease characterised by increased bone fragility and a heightened risk of fractures. It predominantly affects older women, with the prevalence increasing annually [ 13 ]. Osteoporosis arises because of an imbalance between bone resorption and formation. Current treatment options include anti‐resorptive and anabolic drugs; however, many of these medications have side effects and require long‐term use, highlighting the urgent need for novel therapeutic approaches to combat osteoporosis.
Large‐scale randomised clinical trials (RCTs) are a key method for evaluating the efficacy of drug treatments and have increasingly become a valuable tool for integrating human genetics with drug development. RCTs have effectively evaluated the impact of HMGCR and PCSK9 inhibitors on coronary artery disease risk [ 14 , 15 ]. However, RCTs are resource‐intensive and time‐consuming. MR analyses use genetic variants as instrumental variables to determine causal links between exposures and outcomes. Cis‐expression quantitative trait loci (cis‐eQTLs) have been employed to explore the influence of gene expression regulators in diseases such as COVID‐19 and Parkinson's disease [ 16 , 17 ]. Studies have shown that targeting LDL cholesterol may help prevent and manage abdominal aortic aneurysms, benefiting patients [ 18 ].
In this study, we systematically identified VEGFA, PSMA4 and FMO4 as genes significantly associated with osteoporosis risk using MR, and explored their expression patterns using single‐cell RNA sequencing (scRNA‐seq) data from osteoporosis patients. We further investigated their immunological relevance through immune cell infiltration analysis. The integration of genetic, transcriptomic and immunological data provides comprehensive insights into the potential roles of these genes in the pathophysiology of osteoporosis.
FMO4, a member of the flavin‐containing monooxygenase (FMO) family, has not been widely studied in osteoporosis. However, other family members such as FMO3 and FMO5 have been implicated in metabolic and inflammatory regulation—both of which are relevant to bone homeostasis [ 19 ]. Although direct evidence linking FMO4 to bone metabolism is currently limited, previous studies have shown that other FMO family members, such as FMO3 and FMO5, play critical roles in the regulation of metabolism and inflammation—processes that may indirectly influence bone homeostasis. Scott et al. [ 20 ] have suggested that FMO enzymes may influence oxidative stress and immune signalling, affecting osteoblast and osteoclast activity. Our MR analysis identified a significant inverse association between FMO4 expression and osteoporosis risk, suggesting a protective role. Single‐cell analysis revealed enriched expression of FMO4 in vascular smooth muscle cells (VSMCs) and immune cells, indicating a potential role in regulating redox balance and maintaining bone marrow microenvironment stability, thereby influencing bone remodelling indirectly.
PSMA4, which encodes the α4 subunit of the 20S proteasome core, is a key component of the ubiquitin–proteasome system (UPS). This system is involved in protein degradation, immune regulation, and stress responses. Dysregulated proteasome activity has been linked to abnormalities in osteogenesis and bone metabolism [ 21 ]. Clinically, the proteasome inhibitor bortezomib, widely used in multiple myeloma, not only promotes tumour cell apoptosis and suppresses angiogenesis‐related factors but also reduces the invasive effects of myeloma cells on the bone matrix by reversing the hyperactive osteoclast environment, exerting both anti‐myeloma and bone‐protective effects [ 22 ]. Given that proteasomal function influences both osteoclast and osteoblast activity, PSMA4 may serve as a regulator of the osteoclastic/osteogenic balance. In support of this hypothesis, our MR analysis revealed a significant association between increased PSMA4 expression and a higher risk of osteoporosis. These findings highlight PSMA4 as a potential contributor to osteoporosis pathogenesis and a candidate for further mechanistic and therapeutic investigation.
VEGFA, a well‐established member of the vascular endothelial growth factor (VEGF) family, encodes a 27‐kDa dimeric glycoprotein critical to angiogenesis and bone remodelling. In addition to promoting neovascularisation, VEGFA acts as a bone anabolic factor, regulating bone formation and resorption to maintain skeletal homeostasis [ 23 , 24 , 25 ]. In our study, elevated VEGFA expression was significantly associated with an increased risk of osteoporosis, suggesting a potential pathogenic role. This finding is consistent with that of Senel et al. [ 26 ], who reported significantly lower serum VEGF levels in postmenopausal osteoporosis patients compared to healthy controls, supporting VEGFA's importance in bone metabolism. Moreover, Wu et al. [ 1 ] reported that VEGFA expression increases with age. Xu et al. [ 27 ] further demonstrated that oestrogen receptor β (ER‐β) upregulates VEGFA, which subsequently activates the downstream c‐JUN/GADD45G signalling axis, leading to osteoblast apoptosis and accelerated bone loss. These observations are in line with our results and collectively suggest that VEGFA plays a critical role in the pathogenesis of osteoporosis. Taken together, these findings indicate that VEGFA may not only serve as a key regulator of bone homeostasis but also represent a promising therapeutic target for future intervention strategies.
Pathway enrichment analyses indicated that FMO4, PSMA4 and VEGFA may influence osteoporosis progression through several critical signalling pathways, including calcium signalling, Wnt/β‐catenin, PI3K/Akt and Hedgehog signalling. Calcium signalling in osteoclasts is essential for regulating differentiation, bone resorption and transcription. Sustained low intracellular calcium levels can activate nuclear factor of activated T‐cells (NFAT), promoting osteoclastogenesis. The Wnt/β‐catenin pathway plays a central role in maintaining bone homeostasis by promoting osteoblast differentiation and preventing osteoclast apoptosis [ 28 , 29 ]. The PI3K‐Akt‐GSK3β pathway is another important axis in bone disease research. Its activation is closely linked to caspase‐3–mediated apoptosis and can enhance osteoclast activity [ 30 ]. The Hedgehog pathway, in cooperation with bone morphogenetic proteins (BMPs), regulates the differentiation of mesenchymal stem cells into osteoblasts and plays a key role in skeletal development and bone maintenance.
Immune infiltration analysis revealed significant correlations between the three genes and critical immune cells, particularly DCs and macrophages. These immune populations are known to modulate inflammation, bone resorption and remodelling, reinforcing the immunomodulatory relevance of FMO4, PSMA4 and VEGFA in osteoporosis pathophysiology.
A notable strength of this study lies in its focus on druggable genes, which enhances the clinical and translational value of the findings. By integrating genetic, transcriptomic and immunologic data, we identified FMO4, PSMA4 and VEGFA as promising therapeutic targets. Additionally, several small‐molecule inhibitors targeting these genes are currently under development, laying the groundwork for future drug discovery and repurposing.
Nevertheless, this study has limitations. First, the genetic effect sizes used in MR are generally modest, and the estimated causal effects reflect lifelong exposure, which may not directly translate to the outcomes of short‐term therapeutic interventions. Second, the GWAS data used in this study predominantly derive from individuals of European descent, which may limit the generalisability of the findings to other ethnic populations. Future studies should aim to replicate these findings in larger and more diverse cohorts and conduct experimental validation to clarify the biological mechanisms involved.
Conclusions
This study identified key genes linked to osteoporosis, such as FMO4, PSMA4 and VEGFA, and explored their potential molecular mechanisms and regulatory networks. These findings offer valuable insights into the pathological mechanisms underlying osteoporosis and suggest potential candidate genes and pathways for further investigation in drug development and therapeutic research.
Experimental Section
Data Download
Exposure Data: The eQTL data were retrieved from the eQTLGen Consortium database ( https://www.eqtlgen.org ). The eQTLGen Consortium seeks to clarify the genetic structure influencing blood gene expression and enhance comprehension of the genetic foundations of complex traits. The eQTLGen project, now in its second phase, aims to perform extensive genome‐wide meta‐analyses of blood expression data.
Outcome Data: The outcome‐related GWAS data included participants of predominantly European ancestry. The outcome summary data for the training set were obtained from the FinnGen biobank database (finngen_R10_M13_OSTEOPOROSIS). In the FinnGen study, cases were identified using relevant ICD codes. Specifically, the osteoporosis cohort comprised 8017 cases and 391,037 controls. The test set outcome summary data were sourced from the EBI database (GCST90018887). In the EBI study, cases were similarly defined using ICD codes, with the osteoporosis cohort consisting of 18,314 cases and 648,913 controls.
The Gene Expression Omnibus (GEO) database, maintained by the National Center for Biotechnology Information (NCBI), is available at https://www.ncbi.nlm.nih.gov/geo/info/datasets.html . The GEO database is a comprehensive resource for gene expression data. Osteoporosis gene expression data were obtained from the GEO database, specifically from dataset GSE56815 , annotated with the GPL96 platform. This dataset includes 40 control samples and 40 disease samples and was used to explore the bioinformatic mechanisms of key genes in osteoporosis.
The Gene Expression Omnibus (GEO), a public repository maintained by the National Center for Biotechnology Information (NCBI), was utilised to retrieve the single‐cell RNA sequencing (scRNA‐seq) dataset GSE147287 , which comprises profiling data from a single osteoporosis patient.
MR Analysis
The outcome‐related GWAS IDs were extracted from the GWAS summary data ( https://gwas.mrcieu.ac.uk/ ) for causal relationships in eQTLs. SNPs with significant associations ( p < 1e‐5) for each gene across the genome were chosen as potential instrumental variables (IVs). SNPs with R 2 < 0.001 (clumping window size = 10,000 kb) were considered, and weak instrumental variables were excluded by applying an F‐value >10 filter. Four statistical methods were employed to assess causal relationships: Inverse Variance Weighted (IVW), which combines each SNP's Wald estimate using meta‐analysis; MR‐Egger, relying on the Instrument Strength Independent of Direct Effects (InSIDE) assumption; Weighted Median, permitting accurate causal inference even if up to 50% of instrumental variables are invalid; and Weighted Mode, offering greater power than MR‐Egger with reduced bias and lower type I error rates. The Wald ratio method was applied when a single SNP was implicated in the causal relationship. These methods estimated the impact of cis‐acting and certain trans‐acting gene expressions in whole blood on osteoporosis.
Sensitivity Analysis
A leave‐one‐out sensitivity analysis was conducted using MR to evaluate the impact of individual genetic variants on osteoporosis risk. This method sequentially omits each SNP and recalculates the aggregate effect size of the remaining SNPs to detect and eliminate variants that excessively influence the overall estimate. Removing each SNP provides a new point estimate and its 95% confidence interval, allowing assessment of each SNP's unique contribution and the robustness of the results. The summary included both the overall estimate with all SNPs and the estimates after each SNP was excluded. Comparing these estimates allowed us to assess the robustness of our findings by observing the impact of removing individual SNPs on the overall results.
Colocalisation Analysis
Colocalisation analysis was performed using the coloc method, with eQTL summary data and the GWAS data for osteoporosis. A 100‐kilobase region surrounding each index SNP was used to compute posterior probabilities. In the coloc results, H3 represents the posterior probability that gene expression and osteoporosis are linked but have separate causal variants, whereas H4 denotes the probability that they are associated and share a common causal variant. A coloc threshold of SNP.PP.H4 > 0.75 was used to identify shared causal variants.
GSEA
Patients were categorised into high and low expression groups for GSEA on the basis of gene expression levels. Pathway differences between the two groups were analysed using GSEA. The analysis utilised a background gene set sourced from the MsigDB database, offering annotated gene sets for pathways specific to subtypes. Pathway differential expression between subtypes was examined, ranking the most significantly enriched gene sets (adjusted P‐values < 0.05) by their consistency scores. GSEA is commonly used for investigating the biological significance of disease subtypes and their associated pathways.
GSVA
GSVA is an unsupervised, non‐parametric approach for assessing gene set enrichment in transcriptomic datasets. GSVA transforms gene‐level changes into pathway‐level variations by assigning a comprehensive score to each gene set of interest, thereby assessing the biological functions of samples. Gene sets from the Molecular Signatures Database (MSigDB, version 7.0) were analysed using the GSVA algorithm to calculate enrichment scores, facilitating the evaluation of potential biological function variations across samples.
Immune Cell Infiltration Analysis
The single‐sample Gene Set Enrichment Analysis (ssGSEA) method is widely used to evaluate immune cell types in the microenvironment. This method identifies 29 distinct human immune cell phenotypes, such as T cells, B cells, and NK cells. This study utilised ssGSEA to quantify immune cell infiltration from expression profiles, allowing inference of the relative proportions of 29 immune cell types in the samples.
Transcriptional Regulation Analysis of Important Genes
Transcription factor prediction was conducted using the R package ‘RcisTarget.’ RcisTarget's computations rely on motif analysis. The NES for each motif is influenced by the overall count of motifs in the database. Beyond the annotated motifs, additional annotation files were derived using motif similarity and gene sequence analysis. Estimating motif overexpression in gene sets begins with calculating the area under the curve (AUC) for each motif–motif set pair. This is done by calculating the recovery curve on the basis of the ranking of motifs in the gene set. The NES for each motif is calculated on the basis of the AUC distribution of all motifs within the gene set.
miRNA Network Construction
MicroRNAs (miRNAs) are small non‐coding RNAs that control gene expression by facilitating mRNA degradation or blocking mRNA translation. We investigated miRNAs linked to essential genes to determine their role in regulating gene transcription or degradation. Key gene miRNA data were sourced from the miRcode database, and a miRNA‐gene network was visualised with Cytoscape software.
Prediction and Analysis of Drug Interactions With Key Genes
The Drug‐Gene Interaction Database (DGIdb) serves as an extensive resource detailing established and possible gene‐drug interactions. The dataset includes over 14,000 drug‐gene interactions, involving 2600 genes and 6300 drugs, along with 6700 additional genes. By focusing on key genes, we identified potential therapeutic drugs and visualised the gene‐drug interaction network using Cytoscape software.
Quality Control of Single‐Cell Data
Initial preprocessing was conducted using the Seurat package. Cells were filtered on the basis of criteria including total unique molecular identifier (UMI) counts, number of detected genes per cell, and mitochondrial gene expression proportion (defined as the percentage of total gene expression attributable to mitochondrial genes). Aberrant cells exhibiting mitochondrial gene proportions exceeding three median absolute deviations (MADs) from the median were excluded to remove low‐quality or dying cells. DoubletFinder v2.0.4 was employed to identify and eliminate doublet cells, ensuring dataset purity.
Dimensionality Reduction, Clustering and Cell Annotation
Cells underwent global normalisation using the LogNormalize method, scaling gene expression to a target sum of 10,000 transcripts per cell followed by log transformation. Cell cycle scores were computed using CellCycleScoring and highly variable genes (HVGs) were identified via FindVariableFeatures. Confounding factors (e.g., mitochondrial/ribosomal gene expression, cell cycle phase) were regressed out using ScaleData. Principal component analysis (PCA) was applied for linear dimensionality reduction, with principal components (PCs) selected for subsequent analysis. Batch effects were mitigated using the Harmony algorithm, and non‐linear dimensionality reduction was achieved via Uniform Manifold Approximation and Projection (UMAP). Cell types were annotated by integrating marker genes from CellMarker, PanglaoDB, and published literature, supplemented by automated annotation using the SingleR package.
External Validation of the Identified Key Genes
To validate the expression patterns of the hub genes in an independent cohort, we analysed the publicly available microarray dataset GSE35959 (platform: GPL570 ), which contains peripheral blood samples from osteoporosis patients ( n = 5) and healthy controls ( n = 14). The raw data were retrieved from the Gene Expression Omnibus (GEO) database and preprocessed using R software (v4.3.0). Background correction and quantile normalisation were performed using the limma package. Expression values of FMO4, PSMA4 and VEGFA were extracted for comparison. Group differences were assessed using an unpaired Student's t‐test, and results were visualised as box plots generated with the ggplot2 package. Statistical significance was defined as p < 0.05.
Statistical Analysis
To ensure the validity of the MR analysis, we adhered to the three core assumptions: (1) relevance, where instrumental variables (IVs) are strongly associated with the exposure; (2) independence, indicating that IVs are not associated with confounders; and (3) exclusion restriction, meaning that IVs affect the outcome only through the exposure. Violation of the third assumption may indicate the presence of horizontal pleiotropy.
To assess potential horizontal pleiotropy, we performed MR‐Egger regression. A statistically significant non‐zero intercept suggests the presence of directional pleiotropy and potential bias in the causal estimate. In addition, Cochran's Q test was used to evaluate heterogeneity among SNP‐specific causal estimates. A P‐value greater than 0.05 for the Q statistic was considered indicative of no significant heterogeneity.
To further evaluate the robustness of the results and identify any influential single nucleotide polymorphisms (SNPs), a leave‐one‐out sensitivity analysis was conducted. This analysis iteratively removes one SNP at a time to assess its impact on the overall causal estimate.
All statistical analyses were conducted using R software (version 4.3.0). All tests were two‐sided, and a P‐value of less than 0.05 was considered statistically significant.
Author Contributions
Huichao Fu: conceptualization (lead), data curation (supporting), formal analysis (equal), investigation (lead), writing – original draft (lead), writing – review and editing (lead). Yunjiao Wu: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), writing – original draft (equal). Hongfei Lv: methodology (equal), software (equal). Li‐an Qiu: methodology (equal), project administration (equal), resources (equal). Weifeng Hu: resources (equal), software (equal). Hongyun Wu: methodology (equal), project administration (equal). Junda Qian: resources (equal), software (equal). Xiongjie Liang: investigation (equal), methodology (equal), project administration (equal). Xiaoyan Wang: resources (equal), software (equal). Gongping Xu: conceptualization (equal), writing – original draft (supporting), writing – review and editing (supporting).
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgements
This study utilised data from the eQTLGen Consortium ( https://www.eqtlgen.org/ ) and the FinnGen Consortium, which provided summary‐level data on osteoporosis. We acknowledge the investigators of these consortia and studies for generating the data utilised in this analysis.
Fu H., Wu Y., Lv H., et al., “Unravelling Osteoporosis: Key Genes and Potential Therapies,” Journal of Cellular and Molecular Medicine 29, no. 15 (2025): e70759, 10.1111/jcmm.70759.
Huichao Fu and Yunjiao Wu contributed equally as the first author
Contributor Information
Huichao Fu, Email: huichaofu@hrbmu.edu.cn.
Yunjiao Wu, Email: yunjiaowu@hrbmu.edu.cn.
Hongfei Lv, Email: 1796603503@qq.com.
Weifeng Hu, Email: 1638525906@qq.com.
Hongyun Wu, Email: 1094988105@qq.com.
Junda Qian, Email: qianjdsmile@163.com.
Xiongjie Liang, Email: 449004492@qq.com.
Xiaoyan Wang, Email: 5694287@qq.com.
Gongping Xu, Email: 13704846805@163.com.
Data Availability Statement
All data are available in databases. EQTL data were obtained from the eQTLGen consortium database. Osteoporosis GWAS data were obtained from FinnGen biobank and GWAS catalogue database. Osteoporosis microarray data were obtained from GSE56815 .
References
1. Wanionok N. E., Morel G. R., and Fernández J. M., “Osteoporosis and Alzheimer's Disease (Or Alzheimer's Disease and Osteoporosis),” Ageing Research Reviews 99 (2024): 102408.
2. Sietsema D. L., “Fighting the Epidemic: Bone Health and Osteoporosis,” Nursing Clinics of North America 55 (2020): 193–202.
3. Osnes E. K., Lofthus C. M., Meyer H. E., et al., “Consequences of Hip Fracture on Activities of Daily Life and Residential Needs,” Osteoporosis International 15 (2004): 567–574.
4. Kanis J. A., Harvey N. C., Cooper C., Johansson H., Odén A., and McCloskey E. V., “A Systematic Review of Intervention Thresholds Based on FRAX : A Report Prepared for the National Osteoporosis Guideline Group and the International Osteoporosis Foundation,” Archives of Osteoporosis 11 (2016): 25.
5. Anagnostis P., Paschou S. A., Mintziori G., et al., “Drug Holidays From Bisphosphonates and Denosumab in Postmenopausal Osteoporosis: EMAS Position Statement,” Maturitas 101 (2017): 23–30.
6. Lawlor D. A., Harbord R. M., Sterne J. A., Timpson N., and Davey S. G., “Mendelian Randomization: Using Genes as Instruments for Making Causal Inferences in Epidemiology,” Statistics in Medicine 27 (2008): 1133–1163.
7. Davies N. M., Holmes M. V., and Smith D. G., “Reading Mendelian Randomisation Studies: A Guide, Glossary, and Checklist for Clinicians,” BMJ (Clinical Research Ed.) 362 (2018): k601.
8. Ference B. A., Kastelein J., Ray K. K., et al., “Association of Triglyceride‐Lowering LPL Variants and LDL‐C‐Lowering LDLR Variants With Risk of Coronary Heart Disease,” JAMA 321 (2019): 364–373.
9. Klarin D., O'Donnell C. J., and Kathiresan S., “Mendelian Randomization Study of ACLY and Cardiovascular Disease,” New England Journal of Medicine 383 (2020): e50.
10. Zhao J. V., Liu F., Schooling C. M., Li J., Gu D., and Lu X., “Using Genetics to Assess the Association of Commonly Used Antihypertensive Drugs With Diabetes, Glycaemic Traits and Lipids: A Trans‐Ancestry Mendelian Randomisation Study,” Diabetologia 65 (2022): 695–704.
11. Sofat R., Hingorani A. D., Smeeth L., et al., “Separating the Mechanism‐Based and Off‐Target Actions of Cholesteryl Ester Transfer Protein Inhibitors With CETP Gene Polymorphisms,” Circulation 121 (2010): 52–62.
12. Finan C., Gaulton A., Kruger F. A., et al., “The Druggable Genome and Support for Target Identification and Validation in Drug Development,” Science Translational Medicine 9 (2017): eaag1166.
13. Elonheimo H., Lange R., Tolonen H., and Kolossa‐Gehring M., “Environmental Substances Associated With Osteoporosis‐A Scoping Review,” International Journal of Environmental Research and Public Health 18 (2021): 738.
14. Ference B. A., Robinson J. G., Brook R. D., et al., “Variation in PCSK9 and HMGCR and Risk of Cardiovascular Disease and Diabetes,” New England Journal of Medicine 375 (2016): 2144–2153.
15. Cohen J. C., Boerwinkle E., T. H. Mosley, Jr. , and Hobbs H. H., “Sequence Variations in PCSK9, Low LDL, and Protection Against Coronary Heart Disease,” New England Journal of Medicine 354 (2006): 1264–1272.
16. Gaziano L., Giambartolomei C., Pereira A. C., et al., “Actionable Druggable Genome‐Wide Mendelian Randomization Identifies Repurposing Opportunities for COVID‐19,” Nature Medicine 27 (2021): 668–676.
17. Storm C. S., Kia D. A., Almramhi M. M., et al., “Finding Genetically‐Supported Drug Targets for Parkinson's Disease Using Mendelian Randomization of the Druggable Genome,” Nature Communications 12 (2021): 7342.
18. Harrison S. C., Holmes M. V., Burgess S., et al., “Genetic Association of Lipids and Lipid Drug Targets With Abdominal Aortic Aneurysm: A Meta‐Analysis,” JAMA Cardiology 3 (2018): 26–33.
19. Tynes R. E. and Hodgson E., “Catalytic Activity and Substrate Specificity of the Flavin‐Containing Monooxygenase in Microsomal Systems: Characterization of the Hepatic, Pulmonary and Renal Enzymes of the Mouse, Rabbit, and Rat,” Archives of Biochemistry and Biophysics 240 (1985): 77–93.
20. Huang S., Howington M. B., Dobry C. J., Evans C. R., and Leiser S. F., “Flavin‐Containing Monooxygenases Are Conserved Regulators of Stress Resistance and Metabolism,” Frontiers in Cell and Development Biology 9 (2021): 630188.
21. Alexandrova E. M. and Thomsen G. H., “Smurf1 Regulates Neural Patterning and Folding in Xenopus Embryos by Antagonizing the BMP/Smad1 Pathway,” Developmental Biology 299 (2006): 398–410.
22. Dimopoulos M. A., Moreau P., Palumbo A., et al., “Carfilzomib and Dexamethasone Versus Bortezomib and Dexamethasone for Patients With Relapsed or Refractory Multiple Myeloma (ENDEAVOR): A Randomised, Phase 3, Open‐Label, Multicentre Study,” Lancet Oncology 17 (2016): 27–38.
23. Nagy J. A., Dvorak A. M., and Dvorak H. F., “VEGF‐A and the Induction of Pathological Angiogenesis,” Annual Review of Pathology 2 (2007): 251–275.
24. Sağsöz H., Liman N., and Alan E., “Expression of Vascular Endothelial Growth Factor Receptors and Their Ligands in Rat Uterus During the Postpartum Involution Period,” Biotechnic & Histochemistry 90 (2015): 361–374.
25. Wu Y., Cao H., Yang Y., et al., “Effects of Vascular Endothelial Cells on Osteogenic Differentiation of Noncontact Co‐Cultured Periodontal Ligament Stem Cells Under Hypoxia,” Journal of Periodontal Research 48 (2013): 52–65.
26. Senel K., Baykal T., Seferoglu B., et al., “Circulating Vascular Endothelial Growth Factor Concentrations in Patients With Postmenopausal Osteoporosis,” Archives of Medical Science 9 (2013): 709–712.
27. Xu Z., Wang P., Wang Z., et al., “ER‐β Accelerates the Process of Primary Osteoporosis by Promoting VEGFA‐Mediated Apoptosis of Osteoblasts,” Genomics 115 (2023): 110743.
28. Wang C., Chen Q., and Xu H., “Wnt/β‐Catenin Signal Transduction Pathway in Prostate Cancer and Associated Drug Resistance,” Discover Oncology 12 (2021): 40.
29. Han L., Wu J., Wang M., et al., “RNA Modification‐Related Genetic Variants in Genomic Loci Associated With Bone Mineral Density and Fracture,” Genes (Basel) 13 (2022): 1892.
30. Molz S., Dal‐Cim T., Budni J., et al., “Neuroprotective Effect of Guanosine Against Glutamate‐Induced Cell Death in Rat Hippocampal Slices Is Mediated by the Phosphatidylinositol‐3 Kinase/Akt/ Glycogen Synthase Kinase 3β Pathway Activation and Inducible Nitric Oxide Synthase Inhibition,” Journal of Neuroscience Research 89 (2011): 1400–1408.
Associated Data
Data Availability Statement
All data are available in databases. EQTL data were obtained from the eQTLGen consortium database. Osteoporosis GWAS data were obtained from FinnGen biobank and GWAS catalogue database. Osteoporosis microarray data were obtained from GSE56815 .