Introduction
Osteoporosis is characterized by reduced bone strength and increased fracture susceptibility. Bone remodeling involves coordinated osteoblast and osteoclast activity, with contributions from endocrine, metabolic, and immune signals.1–3 Understanding these interactions may help identify molecular candidates for further study.
Experimental studies link microbial metabolites to bone remodeling. Short-chain fatty acids can influence osteoclast differentiation and bone mass, while butyrate can support bone formation through regulatory T-cell-dependent WNT10B signaling.4,5 Studies in germ-free and estrogen-deficient mice further implicate the microbiota and intestinal T-cell trafficking in skeletal homeostasis.6–8 Together, these findings define a gut-bone framework in which microbial products and immune responses can influence both bone formation and resorption.9 These findings motivate investigation of microbiota-associated pathways and their relationship to osteoporosis-associated gene expression.
MR uses genetic instruments to investigate exposure-outcome relationships under assumptions concerning instrument relevance, confounding, and pleiotropy.10,11 Integrating these estimates with gene-expression analyses provides a complementary strategy for prioritizing molecular candidates. Machine-learning approaches provide complementary tools for exploring gene-expression patterns and prioritizing candidate markers.12
The present study integrates multi-omics MR, gene prioritization, public transcriptomic datasets, exploratory machine learning, and RANKL-induced osteoclast differentiation. We aimed to identify candidate genes associated with osteoporosis-related phenotypes and assess whether their protein expression changes during osteoclastogenesis. Combining these approaches enables candidate prioritization at genetic, transcriptomic, and protein levels within a common research framework.
Materials and Methods
Study Design and Multi-Omics Data Acquisition
The analytical workflow is summarized in Figure 1. Osteoporosis GWAS summary statistics were obtained through OpenGWAS under the identifier finn-b-M13_OSTEOPOROSIS, comprising 3,203 cases and 209,575 controls (212,778 participants).13 Gut microbiota data were from the Dutch Microbiome Project (7,738 participants).14 Immune-trait data were from the study by Orrù et al, and metabolite data were from the atlas of 1091 metabolites and 309 metabolite ratios by Chen et al15,16 The corrected MR analysis included 201 microbial, 728 immune, and 1,400 metabolite IVW estimates. Transcriptomic data comprised GSE7429 (20 circulating B-cell samples; 10 high-BMD and 10 low-BMD postmenopausal women), GSE35956 (10 bone-marrow mesenchymal stromal-cell samples; 5 older participants with osteoporosis and 5 middle-aged controls), and GSE56815 (80 monocyte samples; 40 high-BMD and 40 low-BMD women, with 20 premenopausal and 20 postmenopausal participants in each BMD group).17–19 These datasets differ in cell source, age, and phenotype definition; low BMD should not be treated as uniformly equivalent to a clinical osteoporosis diagnosis.Figure 1Overview of the integrative analysis.
The workflow combines public GWAS data, MR and gene-prioritization analyses, three GEO datasets, exploratory expression-based models, and cell experiments. The diagram summarizes the analytical plan and does not establish a causal sequence.A workflow diagram of data integration, multi-omics genetic analysis and machine learning for osteoporosis.Three stacked sections outline an integrative workflow. Top section, Data Acquisition and Integration: Three database icons labeled, GSE7429, GSE35956 and GSE56815, feed by arrows into a funnel icon labeled, Normalization and Batch Correction.
A dashed vertical divider separates this from text reading, GWAS summary statistics for Osteoporosis (finn-b-M13_OSTEOPOROSIS.vcf). Below are three labeled icons: Gut Microbiota, Immune Cells and Metabolites. Middle section, Multi-Omics Genetic Analysis: A heading reads, Mendelian Randomization (MR).
On the left, text reads, SNPs and below it, Instrumental Variables (SNPs). Two arrows point to a box labeled, Exposures and an arrow continues to the word, Osteoporosis. Below, a small label reads, Gene Prioritization, next to a split arrow leading to two items: SNP-to-Gene Mapping and SMR (Summary-data-based MR).
These converge toward a two-circle overlap graphic labeled, Key Intersection Genes. Bottom section, Candidate Evaluation and Machine Learning: On the left, a cluster labeled, Immune Signature Patterns, includes labeled cell types: macrophage, T-cell, dendrite, monocyte and mast cell. Arrows lead to four connected circles labeled, GLM, SVM, XGB and RF, with the caption below, Machine Learning Evaluation.
An arrow points to a small chart-like graphic labeled, Exploratory Five-Gene Nomogram.
Two-Sample Mendelian Randomization and Bayesian Colocalization Analysis
Two-sample MR was used to examine associations between systemic exposures and osteoporosis.20 Genetic instruments were selected at exposure-association P < 1 × 10−5 and F > 10.21 The original specification used European-reference LD clumping at r2 < 0.001 within 10,000 kb. During revision, harmonization and inverse-variance-weighted (IVW) estimation were repeated from the exposure and outcome instrument tables, retaining outcome-associated variants rather than excluding them on the basis of outcome P values. Alleles were aligned, and palindromic variants with ambiguous allele frequencies were excluded.
Multiplicative random-effects IVW estimates were used, with residual variance not permitted to reduce standard errors below their fixed-effect values. Benjamini-Hochberg adjustment was applied to all IVW P values within each exposure domain (201, 728, and 1,400 tests); an additional adjustment pooled all 2,329 tests.22 Adjusted P < 0.05 defined the multiplicity threshold. MR-Egger, weighted-median estimation with 10,000 parametric bootstrap replicates, and leave-one-SNP-out IVW analyses were performed for the association meeting this threshold.23,24 Colocalization assessed support for a shared association signal under the specified model.25 Posterior probabilities from the complete colocalization outputs were evaluated; PP.H4 was compared with 0.7 and, descriptively, 0.5.
Summary-Data-Based MR (SMR) and Transcriptomic Harmonization
SMR was used for candidate prioritization with cis-eQTL data.26 The analysis command used SMR version 1.3.1, eQTLGen cis-eQTL summary data,27 a 1000 Genomes European LD reference,28 and a minor allele frequency threshold of 0.01. Rows with missing values were removed before retaining nominal P_SMR < 0.05 and P_HEIDI > 0.05. A post hoc Benjamini-Hochberg assessment was also applied across all 15,666 available SMR P values.
The three normalized GEO matrices were merged over 13,236 shared genes and 110 samples (55 per analysis group). ComBat used study as the batch variable without a phenotype covariate.29 Expression distributions and PCA were inspected before and after adjustment. SNP-mapped genes were intersected with the nominal SMR gene set.
Of the 43 intersection genes, 33 were present in the merged expression matrix. These genes were screened using two-sided Wilcoxon rank-sum tests with an unadjusted P < 0.05; no fold-change threshold was applied. This nominal screen was used for exploratory candidate prioritization.
Post hoc Benjamini-Hochberg-adjusted P values across the 33 tests are also reported. The resulting nine-gene matrix was used for expression, co-expression, and classifier analyses.
Exploratory Immune Signature and Correlation Analysis
CIBERSORT was run with the LM22 reference, 1,000 permutations and quantile normalization enabled.30 The implementation transformed the input using 2 raised to the expression value when its maximum was below 50. Samples were retained at CIBERSORT P < 0.05; the exported fraction matrix contained 110 samples (55 per analysis group). Group comparisons used two-sided Wilcoxon rank-sum tests, and associations with the nine candidate genes used Spearman correlations.
These P values describe exploratory comparisons; post hoc Benjamini-Hochberg adjustments across the 22 group comparisons and, separately, across the estimable gene-signature correlations (198 pairs assessed; constant signatures were not testable) are provided. The inputs represent heterogeneous purified cell populations and mesenchymal stromal cells, so the estimates were interpreted as exploratory LM22 signatures rather than measurements of leukocyte infiltration in bone.
Exploratory Classifier and Nomogram Construction
Nine candidate genes were evaluated in 110 samples using the caret framework.31 With a random seed of 123456, createDataPartition specified a stratified split with training fraction 0.7 (78 training samples, 39 per group; 32 testing samples, 16 per group). Random Forest, radial-kernel support vector machine, xgbDART, and binomial generalized linear models were fitted using a repeated cross-validation control with five folds. Group status (Type) was the response in the model formula.
ROC curves were calculated with pROC,32 and residual summaries were assessed on the test partition. DALEX permutation importance on that partition was subsequently used to select five RF candidates.33 The RF-selected five genes were subsequently used in a separate multivariable logistic regression model fitted with rms: lrm to all 110 samples. A nomogram was generated from this logistic model.34 Calibration used 1,000 bootstrap resamples, and decision-curve analysis was exploratory.35 These procedures did not repeat the preceding gene-selection process within each bootstrap resample and did not constitute external validation.
Additional candidate-gene expression comparisons in GSE2208 and GSE7158 are described in the Supplementary Methods and Results.36,37 External prediction analyses used reconstructed nine-gene RF, five-gene logistic-regression and five-gene RF models, fitted exclusively to the development data. Cohort eligibility, sample overlap, predictor coverage and all direct-transfer and frozen-reference sensitivity results are reported in Supplementary Tables S10–S12 and Supplementary Figures S1, S2. A supplementary reanalysis started from the supplied gene-level matrices before our additional normalization.
Reference quantile normalization,38 nominal screening of the 33 measurable genetically prioritized genes, and RF tuning were fitted within training folds only.39 Three repetitions of fivefold outer validation with threefold inner tuning and three study-held-out analyses were performed without joint ComBat.40 This new Python implementation used scikit-learn41 to evaluate a screening-and-classification pipeline rather than the historical fitted RF or a fixed five-gene panel. Upstream GEO processing, including RMA in GSE56815 and GSE7429, was retained. Full methods, partitions and results are provided in Supplementary Tables S13–S15 and Figure S3.
Cohort-specific descriptive associations for the five original candidates are reported in Supplementary Table S16 and Figure S4. A separate post hoc sensitivity analysis fixed the original five genes without further feature selection. Ridge logistic regression, shrinkage linear discriminant analysis, radial-kernel SVM and RF were compared across three expression representations and three sample populations using the same three repetitions of stratified fivefold cross-validation.
Scaling and reference normalization were fitted within training folds. These analyses were conditional on the prior five-gene selection and did not independently validate that selection. All 36 specifications are reported in Supplementary Tables S17 and S18.
We subsequently assessed ridge logistic-regression models using the same fixed five genes across seven external datasets, including the newly identified GSE56116 blood and GSE230665 femur datasets. The model used L2 regularization (C = 1), with feature scaling and coefficients fitted only in development data. Six specifications combined two development populations (all 110 samples or the 80 monocyte samples) with source log-scale, common-gene rank, or training-reference quantile representations.
The primary specification used all 110 samples and source log-scale values. Rank and quantile variants used 9,249 genes shared across platforms for preprocessing; only the original five genes entered the classifier. External outcomes were not used for fitting or preprocessing, and the prespecified score direction was retained.
These new variants do not constitute direct validation of the earlier 13,236-gene-rank internal model or the historical RF. Cohort screening, complete performance estimates and individual predictions are supplied in Supplementary Tables S21–S23 and Figure S5.
Cell Culture and in vitro Osteoclast Differentiation
The murine macrophage cell line RAW264.7 was obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM; Servicebio, Wuhan, China) supplemented with 10% fetal bovine serum (FBS; Excell Bio, Taicang, China) and 1% penicillin-streptomycin at 37 °C in a humidified atmosphere containing 5% CO2. For osteoclast differentiation, RAW264.7 cells were seeded into 6-well plates.
After 24 hours of adherence, the culture medium was replaced. The cells were divided into two groups: the Control group, which was maintained in standard complete DMEM, and the RANKL+ group, which was cultured in complete DMEM supplemented with 50 ng/mL recombinant murine RANKL (Novoprotein, Suzhou, China). The induction medium was refreshed every two days.
After 5 days of continuous stimulation, cells were processed for downstream analyses. Three independent experiments were performed, with three parallel wells per group in each experiment. The independent experiment was the statistical unit for the protein and TRAP analyses.
Tartrate-Resistant Acid Phosphatase (TRAP) Staining
To confirm the successful induction of mature osteoclasts, cells were washed with PBS and fixed with 4% paraformaldehyde for 15 minutes at room temperature. Subsequently, the cells were stained using a TRAP leukocyte diagnostic kit (Solarbio, Beijing, China) according to the manufacturer’s instructions. TRAP-positive (purplish-red) multinucleated cells containing three or more nuclei were classified as mature osteoclasts.
Images were captured using an inverted optical microscope (Olympus, Tokyo, Japan), and the number of osteoclasts was quantified in five randomly selected, non-overlapping fields per well. Each of the three TRAP values per group represents one independent experiment, rather than an individual field or a technical replicate. Counts from the sampled fields and parallel wells were averaged to obtain one TRAP value for each independent experiment.
Protein Extraction and Western Blot Analysis
Total proteins were extracted from both the Control and RANKL-induced RAW264.7 cells using RIPA lysis buffer (Beyotime, Shanghai, China) supplemented with protease and phosphatase inhibitor cocktails (Servicebio, Wuhan, China). Protein concentrations were determined using a BCA Protein Assay Kit (Beyotime). Equal amounts of protein (30 μg/lane) were separated by 10% SDS-PAGE and electro-transferred onto PVDF membranes (Millipore, Billerica, MA, USA).
After blocking with 5% non-fat milk for 1 hour at room temperature, the membranes were incubated overnight at 4 °C with specific primary antibodies. The primary antibodies used were: anti-CCR1 (1:1000, Proteintech, Wuhan, China), anti-FADS2 (1:1000, Proteintech, Wuhan, China), anti-FBXL14 (1:1000, Affinity Biosciences, Jiangsu, China), anti-TRMT2A (1:1000, ABclonal, Wuhan, China), anti-TTLL12 (1:1000, Bioss, Beijing, China), and anti-GAPDH (1:5000, Proteintech, Wuhan, China) as the internal loading control. Following rigorous washing, membranes were incubated with HRP-conjugated secondary antibodies (1:3000, Proteintech, Wuhan, China) for 1 hour at room temperature.
Protein bands were visualized utilizing an enhanced chemiluminescence (ECL) detection system (Vazyme, Nanjing, China) and quantified by densitometry using ImageJ software (NIH, Bethesda, MD, USA).42 Protein-expression quantification included three independent experiments per group. Control and RANKL-treated groups were compared using unpaired, two-tailed Student’s t-tests with a pooled variance estimate. Descriptive results are summarized as mean ± SD.
Primary experimental P values are unadjusted, with P < 0.05 used as the nominal significance threshold. The sample size for each comparison was three independent experiments per group. Post hoc Benjamini-Hochberg adjustment covered the five protein outcomes.
Two-sided Welch tests with a separate five-protein adjustment assessed sensitivity to unequal variances; TRAP was treated as a separate differentiation readout.
Results
Transcriptomic Data Harmonization
The merged expression matrix contained 13,236 genes and 110 samples, with 55 samples in each analysis group. Figure 2 presents expression distributions and PCA before and after ComBat adjustment, documenting the preprocessing of the three development datasets.Figure 2Expression distributions and PCA before and after ComBat adjustment. (A and B) Sample expression distributions before and after adjustment.
(C and D) Corresponding PCA projections. The input comprises 110 samples: GSE35956, n = 10; GSE56815, n = 80; and GSE7429, n = 20. Ellipses describe the 95% bivariate normal dispersion region for each study.
These panels are descriptive; the before/after comparison does not establish removal of biological differences between the studies.A composite of two boxplot charts and two scatter plots comparing data before and after batch correction.A) Boxplot chart titled, Before batch correction. The x axis label is, Sample, with unit not shown. The y axis label is, Expression, with unit not shown.
The y axis shows tick labels at 5, 10 and 15, giving a visible labeled range from 5 to 15. Many adjacent boxplots are shown across the sample axis. A legend at right is titled, Project, listing GSE35956, GSE56815 and GSE7429.
The boxplots are arranged in three contiguous sample blocks that match the three legend entries. The first block spans most samples and has boxes centered between the 5 and 10 ticks. The second block is shorter and has boxes centered slightly above the 5 tick.
The third block is the shortest and has boxes centered near the 10 tick, with whiskers extending above the 15 tick. B) Boxplot chart titled, After batch correction. The x axis label is, Sample, with unit not shown.
The y axis label is, Expression, with unit not shown. The y axis shows tick labels at 3, 6, 9 and 12, giving a visible labeled range from 3 to 12. Many adjacent boxplots are shown across the sample axis.
A legend at right is titled, Project, listing GSE35956, GSE56815 and GSE7429. The three contiguous sample blocks align more closely in median height than in A), with most boxes centered between the 6 and 9 ticks and whiskers extending down toward 3 and up toward 12. C) Scatter plot titled, Before batch correction.
The x axis label is, PC1, with unit not shown. The y axis label is, PC2, with unit not shown. The x axis shows tick labels at negative 300, negative 200, negative 100 and 0, giving a visible labeled range from negative 300 to 0.
The y axis shows tick labels at negative 30, 0, 30, 60 and 90, giving a visible labeled range from negative 30 to 90. A legend at right is titled, Type, listing GSE35956, GSE56815 and GSE7429, each with a distinct point marker. Three separated point clusters appear: GSE35956 points cluster near PC1 around negative 300 and PC2 around 0 to about 10; GSE56815 points cluster near PC1 around 0 and PC2 around negative 30 to about negative 10; GSE7429 points cluster near PC1 around 0 and PC2 around about 60 to near 90.
D) Scatter plot titled, After batch correction. The x axis label is, PC1, with unit not shown. The y axis label is, PC2, with unit not shown.
The x axis shows tick labels at negative 30, negative 20, negative 10, 0, 10, 20 and 30, giving a visible labeled range from negative 30 to 30. The y axis shows tick labels at negative 10, 0, 10 and 20, giving a visible labeled range from negative 10 to 20. A legend at right is titled, Type, listing GSE35956, GSE56815 and GSE7429, each with a distinct point marker.
Points from all three types overlap in a single central cloud, concentrated near PC1 from about negative 20 to 20 and PC2 from about negative 10 to 15, with a few points above the 20 tick on PC2. Three ellipses are drawn around the overlapping point cloud, one ellipse per type, centered near PC1 around 0 and PC2 around 0. The four graphs pair distribution views in A) and B) with projection views in C) and D) for the same before and after batch correction conditions.
Multi-Omics MR Associations
After retaining instruments regardless of outcome P value, IVW analyses comprised 201 microbial, 728 immune, and 1,400 metabolite tests. Nominal P < 0.05 occurred for 8, 54, and 64 traits, respectively. One metabolite association survived Benjamini-Hochberg adjustment: the metabolite exposure indexed by GCST90200621 (26 SNPs; OR 1.296, 95% CI 1.155–1.454; P = 9.64 × 10−6; within-domain adjusted P = 0.0135; pooled adjusted P = 0.0224).
No microbial or immune association passed adjustment. For this metabolite association, the MR-Egger and weighted-median estimates had the same direction (P = 0.00594 and 0.0124, respectively); the Egger intercept P was 0.427. Leave-one-SNP-out IVW odds ratios ranged from 1.267 to 1.322.
Figure 3 displays corrected estimates for the originally prioritized exposures, including the corrected metabolite association; Supplementary Table S1 reports every corrected IVW estimate. The original nominal shortlist comprised 5 microbial, 14 immune, and 24 metabolite traits and was used for exploratory downstream candidate mapping. Sensitivity estimates and leave-one-SNP-out results for this metabolite association are provided in Supplementary Table S7.Figure 3Corrected IVW estimates for the originally prioritized exposures.
(A) Microbial traits; (B) immune traits; (C) metabolites. Points and horizontal bars indicate odds ratios and 95% confidence intervals after retaining outcome-associated SNPs during harmonization. Red and cyan squares indicate odds ratios above and below 1, respectively; bold nominal P values are below 0.05.
Multiplicity-adjusted results are reported in Supplementary Table S1. The displayed shortlist was defined by the original nominal screen, rather than reselected from the corrected P values. The corrected metabolite association meets the corrected multiplicity threshold.
Supplementary Table S1 provides the SNP counts, nominal P values, and within-domain and pooled adjusted P values for all 2329 tested traits.Three forest plots of inverse variance weighted odds ratios for microbial, immune and metabolite traits.A) Forest plot table of microbial traits. Columns are labeled, exposure, nsnp, method, pval and OR with 95 percent CI. The method shown for all rows is Inverse variance weighted.
The x-axis label is not shown; the plotted measure is odds ratio with 95 percent confidence interval. Five exposures are listed with point estimates and horizontal confidence intervals: o underscore Lactobacillales f underscore Lactobacillaceae, nsnp 12, pval 0.015, odds ratio 1.124 with 95 percent confidence interval 1.023 to 1.234; g underscore Lactobacillaceae g underscore Lactobacillus, nsnp 11, pval 0.017, odds ratio 1.124 with 95 percent confidence interval 1.021 to 1.237; g underscore Bacteroides s underscore Bacteroides massiliensis, nsnp 6, pval 0.008, odds ratio 1.298 with 95 percent confidence interval 1.069 to 1.576; g underscore Bacteroides s underscore Bacteroides stercoris, nsnp 7, pval 0.014, odds ratio 1.300 with 95 percent confidence interval 1.055 to 1.601; g underscore Coprococcus s underscore Coprococcus comes, nsnp 9, pval 0.015, odds ratio 0.804 with 95 percent confidence interval 0.674 to 0.959. B) Forest plot table of immune traits with the same columns, measure and method as A.
Fourteen exposures are listed: IgD minus CD38dim AC, nsnp 14, pval 0.201, odds ratio 1.103 with 95 percent confidence interval 0.949 to 1.281; CD62L minus DC AC, nsnp 13, pval 0.007, odds ratio 1.151 with 95 percent confidence interval 1.040 to 1.274; CD62L minus myeloid DC AC, nsnp 15, pval 0.002, odds ratio 1.127 with 95 percent confidence interval 1.043 to 1.218; EM DN CD4 minus CD8 minus T cell, nsnp 20, pval 0.047, odds ratio 1.098 with 95 percent confidence interval 1.001 to 1.204; T slash B, nsnp 24, pval 0.077, odds ratio 0.936 with 95 percent confidence interval 0.869 to 1.007; CD4 plus CD8dim leukocyte, nsnp 14, pval 0.010, odds ratio 0.873 with 95 percent confidence interval 0.787 to 0.969; CD28 plus CD45RA minus CD8dim AC, nsnp 15, pval 0.00028, odds ratio 0.840 with 95 percent confidence interval 0.765 to 0.923; CD28 plus CD127 minus CD25 plus CD8br s T cell, nsnp 16, pval 0.190, odds ratio 1.101 with 95 percent confidence interval 0.953 to 1.271; BAFF R on IgD minus CD24 minus, nsnp 19, pval 0.034, odds ratio 0.937 with 95 percent confidence interval 0.882 to 0.995; BAFF R on IgD minus CD24 minus, nsnp 21, pval 0.042, odds ratio 0.942 with 95 percent confidence interval 0.890 to 0.998; CD38 on naive mature B cell, nsnp 21, pval 0.097, odds ratio 1.066 with 95 percent confidence interval 0.988 to 1.150; HVEN on CM CD8br, nsnp 16, pval 0.007, odds ratio 0.913 with 95 percent confidence interval 0.855 to 0.975; FSCA on monocyte, nsnp 16, pval 0.324, odds ratio 1.061 with 95 percent confidence interval 0.944 to 1.192; HLA DR on HLA DR plus CD8br, nsnp 18, pval 0.094, odds ratio 0.917 with 95 percent confidence interval 0.829 to 1.015. C) Forest plot table of metabolite traits with the same columns, measure and method as A.
Prioritization of Candidate Genes
The SMR results and chromosomal distributions are displayed in Figure 4A and B. The complete colocalization tables contained 58 microbiota-associated, 403 immune-trait-associated, and 823 metabolite-associated records. Maximum PP.H4 values were 0.0453, 0.1816, and 0.2438, respectively; no record exceeded either 0.5 or 0.7 (Figure 4C and Supplementary Table S4).
Separately, SNP-to-gene mapping identified 45 microbial, 338 immune, and 742 metabolite-associated gene symbols, yielding a union of 1,089 genes. The nominally filtered SMR table contained 774 unique genes. None of the 15,666 SMR tests survived the additional Benjamini-Hochberg adjustment (minimum adjusted P = 0.1233).
Their intersection comprised 43 genes, leaving 1,046 genes unique to SNP mapping and 731 unique to SMR (Figure 4D and Supplementary Table S3). These intersections were used for candidate prioritization. The complete SMR results and adjusted P values are provided in Supplementary Table S6.Figure 4Genetic and transcriptomic candidate prioritization.
(A and B) SMR results and chromosomal distributions, highlighting nominally selected genes (P_SMR < 0.05 and P_HEIDI > 0.05); no SMR test passed the post hoc false-discovery-rate threshold. (C) Schematic Venn diagram of genes supported by colocalization at PP.H4 > 0.7: all regions contain zero genes. Circle areas are illustrative.
The complete microbiota, immune-trait and metabolite outputs contain 58, 403 and 823 records, respectively; their maximum PP.H4 values are 0.0453, 0.1816 and 0.2438. No record exceeds 0.5 or 0.7. These record counts are not participant sample sizes or independent loci.
(D) Intersection of 1,089 unique SNP-mapped genes and 774 SMR genes, with 43 shared genes, 1,046 SNP-mapped-only genes, and 731 SMR-only genes (Supplementary Table S3). (E) Expression distributions of the nine nominally selected genes; Wilcoxon P values are in Supplementary Table S5. (F) Row-standardized expression heatmap.
(G) SMR probe coordinates for the nine genes. (H) Spearman co-expression coefficients. Expression analyses use the merged 110-sample matrix.
PP.H4 is a posterior probability, not a P value. Complete colocalization results are provided in Supplementary Table S4. In (E) blue denotes controls and red denotes the comparison group.
In (E) and (F), * indicates nominal Wilcoxon P < 0.05 and ** indicates nominal Wilcoxon P < 0.01; these symbols do not denote multiplicity-adjusted significance.Different plots summarizing SMR, gene overlaps, expression distributions, probe locations and co expression.A) Scatter plot of SMR results across chromosomes. The x axis label is Chr, with chromosome groups labeled 1 through 22. The y axis label is minus log subscript 10 end subscript, P value.
The y axis tick labels are 0, 1, 2, 3, 4, 5, 6. Points form vertical bands within each chromosome group. A dashed horizontal reference line crosses the plot at a y value between 1 and 2.
Several points are annotated with gene names; many annotations are too small to read clearly. B) Circular chromosomal distribution plot. The x axis label and unit are not applicable.
The y axis label and unit are not applicable. The outer ring is segmented and labeled with chromosome numbers 1 through 22 around the circle. Inner rings contain many radial marks and points; individual numeric values are not labeled.
C) Three circle Venn diagram labeled Gut microbiota, Immune cells and Metabolites, with text PP. H4 greater than 0.7 above. Each of the seven regions shows the value 0.
D) Two circle Venn diagram labeled SNP and SMR. The left only region shows 1046 with 57.5 percent. The overlap shows 43 with 2.4 percent.
The right only region shows 731 with 40.2 percent. E) Box and whisker plots of expression for nine genes, with a legend labeled Type containing Control and Osteoporosis. The x axis label is Gene, with gene labels CCR3, ABCC2, FADS2, TRIM72, TTLL12, SMARCA4, FBXL14, CCR1 and FADS1.
The y axis label is Expression, with tick labels 0, 2, 4, 6, 8, 10. Each gene has two boxplots, one per group, with scattered points around them. F) Heatmap with hierarchical clustering at left and a top annotation bar labeled Type.
A legend labeled Type lists Control and Osteoporosis. A vertical scale bar at right shows values 3, 2, 1, 0, negative 1, negative 2, negative 3. Row labels list CCR3, ABCC2, FADS2, TRIM72, TTLL12, SMARCA4, FBXL14, CCR1 and FADS1.
G) Circular plot of SMR probe coordinates. The x axis label and unit are not applicable. The y axis label and unit are not applicable.
Chromosome segments labeled 1 through 22 form an outer ring. Several inner tracks show short arc segments and radial bars at specific genomic positions; numeric coordinates are not labeled. H) Circular chord diagram of Spearman co expression coefficients.
The x axis label and unit are not applicable. The y axis label and unit are not applicable. A vertical scale at right is labeled with 1 at the top, 0 in the middle and negative 1 at the bottom.
Multiple chords connect segments around the circle; individual pair labels are not shown.
The transcriptomic analyses displayed nine candidate genes: CCR1, CCR3, TRMT2A, TTLL12, SMARCA4, FADS1, FADS2, FBXL14, and ABCC2 (Figure 4E–H). The expression plots, chromosomal positions, and co-expression patterns characterize these candidates across the merged transcriptomic data. The nine genes met the nominal Wilcoxon threshold of P < 0.05.
The full 33-gene screening results, including post hoc adjusted P values, are reported in Supplementary Table S5; the nine genes are retained as exploratory candidates rather than multiplicity-confirmed transcriptomic discoveries.
Exploratory Immune Signature Patterns
The CIBERSORT output and gene-signature correlations are shown in Figure 5. These exploratory LM22 patterns describe variation across the selected cell populations and provide context for the candidate-gene expression findings. Fractions, group comparisons and gene-signature correlations are provided in Supplementary Table S8.Figure 5Exploratory LM22 signature estimates and correlations.
(A) Fractions in 110 samples retained by the CIBERSORT significance filter (55 per analysis group). (B) Spearman correlations between nine candidate genes and the 20 estimable LM22 signatures. (C) Signature distributions by analysis group.
Wilcoxon rank-sum comparisons and correlations are exploratory; nominal and Benjamini-Hochberg-adjusted P values are reported in Supplementary Table S8. The color scale in (B) indicates correlation coefficients. The heterogeneous purified-cell inputs do not permit interpretation as measured immune-cell infiltration in bone.
The two constant signatures (T cells follicular helper and Mast cells activated) have undefined correlations and are omitted from (B) and from multiplicity adjustment. In (C), blue denotes controls and red denotes the comparison group.Different plots showing LM22 immune signatures, gene correlations and groupwise signature distributions.A) A stacked bar chart of immune cell type composition by sample for two groups. The vertical axis label is Relative Percent, unit is percent, ranging from 0 percent to 100 percent.
The horizontal axis shows individual samples; tick labels are not legible. A group strip under the bars labels Control and Osteoporosis. Each vertical bar is subdivided into stacked segments corresponding to the legend entries: B cells naive, B cells memory, Plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, T cells gamma delta, NK cells resting, NK cells activated, Monocytes, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Mast cells activated, Eosinophils, Neutrophils.
Across nearly all samples, one dominant stacked component occupies most of each bar, with the remaining cell types forming thinner bands near the bottom and top. Individual numeric segment values are not printed. B) A heatmap of correlations between candidate genes and LM22 signatures.
The horizontal axis lists immune signatures: B cells naive, B cells memory, Plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells regulatory, T cells gamma delta, NK cells resting, NK cells activated, Monocytes, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Eosinophils, Neutrophils and Mast cells activated. The vertical axis lists genes: TTLL12, TRMT2A, SMARCA4, FBXL14, FADS2, FADS1, CCR3, CCR1 and ABCC2. A legend at the right is labeled Correlation with a numeric scale marked at 0.2, 0.0, negative 0.2 and negative 0.4.
A significance key shows asterisks with thresholds: three asterisks p less than 0.001, two asterisks p less than 0.01, one asterisk p less than 0.05. Several heatmap cells contain one, two, or three asterisks; exact per cell correlation numbers are not printed. C) A grouped distribution plot with points and boxplots comparing two groups across immune signatures.
The vertical axis label is Fraction, unit not shown, ranging from 0.0 to 0.8 with labeled ticks at 0.0, 0.2, 0.4, 0.6 and 0.8. The horizontal axis shows categories: B cells naive, B cells memory, Plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, T cells gamma delta, NK cells resting, NK cells activated, Monocytes, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Mast cells activated, Eosinophils, Neutrophils. A legend labeled Type lists Control and Osteoporosis.
For most categories, points cluster close to 0.0 with small spreads. Monocytes show the largest fractions, with boxplots centered around roughly 0.6 to 0.7 and points extending from roughly 0.5 up to roughly 0.8. A single high point near 0.8 appears above the main cluster at the Monocytes category.
Exploratory Classifier Performance
The original test-partition analysis yielded AUCs of 0.805 for RF, 0.793 for SVM, 0.762 for XGB, and 0.785 for GLM (Figure 6B). RF had the highest displayed point estimate in this nine-gene comparison. Its importance output prioritized CCR1, TRMT2A, FADS2, FBXL14, and TTLL12.
Residual summaries and importance output are presented in Figure 6A, C, and D. These comparisons summarize internal model performance in the development data.Figure 6Exploratory comparison of four classifiers. (A) Residual boxplots.
(B) Original internal test ROC curves and reported AUCs: RF, 0.805; SVM, 0.793; XGB, 0.762; GLM, 0.785. (C) Feature-importance output. (D) Residual distributions.
The input comprises nine candidate genes and 110 samples; ROC and residual assessments used the 32-sample test partition (16 per group), following allocation of 78 samples to training (39 per group) with createDataPartition at p = 0.7. No statistical superiority of one AUC over another is established by these point estimates. The Type label in the source importance visualization denotes the group response, which was not an input predictor in the model formulas.
Pre-split screening and subsequent use of test-partition importance limit this estimate; the revised pipeline assessments are reported in Supplementary Figure S3.Different plots comparing four classifiers using residuals, receiver operating curves and feature importance.A) Boxplots of absolute residual. Text above plot reads, Boxplots of absolute residual and, Red dot stands for root mean square of residuals. Legend text reads, Model, GLM, SVM, RF, XGB.
The horizontal axis label is, absolute residual, with no unit shown, ranging from 0.00 to 1.00, with tick labels 0.00, 0.25, 0.50, 0.75, 1.00. The vertical axis lists categories GLM, SVM, RF, XGB, with no axis label and no unit shown. Each category has a horizontal boxplot and a single dot.
GLM box spans approximately 0.15 to 0.55 with the dot near 0.45. SVM box spans approximately 0.20 to 0.60 with the dot near 0.45. RF box spans approximately 0.25 to 0.50 with the dot near 0.40.
XGB box spans approximately 0.10 to 0.60 with the dot near 0.45. B) Receiver operating characteristic curve plot with four step curves and a diagonal reference line. The horizontal axis label is, 1 minus Specificity, with no unit shown, ranging from 0.0 to 1.0.
The vertical axis label is, Sensitivity, with no unit shown, ranging from 0.0 to 1.0. A legend inside the plot lists, RF: 0.805, SVM: 0.793, XGB: 0.762, GLM: 0.785. The four curves rise in steps from near the origin toward the top right corner, with multiple vertical and horizontal segments; the diagonal reference line runs from bottom left to top right.
C) Feature importance horizontal bar charts grouped by model. Title text reads, Feature Importance and below it, created for the GLM, RF, SVM, XGB model. Four stacked sections are labeled GLM, RF, SVM, XGB.
The vertical axis in each section lists the same feature names: FADS2, TTLL12, FBXL14, ABCC2, TRMT2A, CCR3, FADS1, CCR1, SMARCA4. A label, Type, appears near the lower left of each section. The horizontal axis label at the bottom reads, Root mean square error RMSE loss after permutations, with no unit shown.
The horizontal axis shows tick labels 0.42, 0.44, 0.46, 0.48, 0.50. Bars extend rightward from a baseline near 0.42. In GLM, the longest bar is for FADS2 reaching close to 0.49, with other bars of varying shorter lengths.
In RF, the longest bar is for FADS2 reaching close to 0.48, with several mid length bars and multiple short bars near the baseline. In SVM, the longest bar is for FADS2 reaching around 0.46, with TTLL12 and FBXL14 also extending to the mid 0.45 range and several shorter bars near 0.42 to 0.44. In XGB, the longest bar is for FADS2 reaching close to 0.50, with TTLL12 reaching around 0.48 and additional bars extending to roughly 0.44 to 0.46.
D) Reverse cumulative distribution of absolute residual. Title text reads, Reverse cumulative distribution of absolute residual. Legend text reads, Model, GLM, SVM, RF, XGB.
The horizontal axis label is, absolute residual, with no unit shown, ranging from 0.00 to 1.00. The vertical axis shows percentages from 0 percent to 100 percent, with no axis label and no unit beyond percent. Four descending step lines start near 100 percent at low absolute residual values and decrease toward 0 percent as absolute residual approaches 1.00, with the lines closely overlapping across most of the range.
The five prioritized genes were CCR1, TRMT2A, FADS2, FBXL14, and TTLL12. Their individual AUCs ranged from 0.609 to 0.651 (Figure 7A), indicating modest discrimination in the analyzed data. Calibration and decision curves and the separate five-gene logistic-regression nomogram are shown in Figure 7B–D.
Figure 8 summarizes the working biological hypothesis linking the study findings. The external expression comparisons did not reach nominal significance (Supplementary Table S9 and Supplementary Methods and Results). In the primary external monocyte cohort, GSE56814 (n = 73), the reconstructed nine-gene RF yielded an AUC of 0.503 (95% CI 0.359–0.642).
None of the 26 evaluable external model comparisons reached nominal P < 0.05 (Supplementary Tables S10–S12; Supplementary Figures S1 and S2). In the supplementary reanalysis from pre-normalization gene matrices, out-of-fold AUCs were 0.558, 0.521 and 0.505 across three nested-validation repetitions (mean 0.528). Study-held-out AUCs were 0.480 for GSE35956, 0.504 for GSE56815 and 0.510 for GSE7429 (Supplementary Table S13 and Figure S3).
All fold-level features and predictions are supplied. In the original development cohorts, three of 15 five-gene expression comparisons were nominally significant, but none met the Benjamini-Hochberg threshold (minimum q = 0.0756; Supplementary Table S16 and Figure S4). With the original five genes fixed, ridge logistic regression using within-sample expression percentiles yielded internal AUCs of 0.744–0.766 (mean 0.756) across the three repetitions in 110 samples, with Brier scores of 0.202–0.206 compared with 0.250 for a prevalence-only prediction.
Training-reference quantile normalization yielded ridge-model AUCs of 0.744–0.761 (mean 0.753). These findings support moderate internal discrimination conditional on the selected panel; the complete post hoc comparison is reported in Supplementary Tables S17 and S18.Figure 7Exploratory five-gene expression model. (A) Single-gene ROC curves: CCR1, AUC 0.615; TRMT2A, 0.609; FADS2, 0.623; FBXL14, 0.651; TTLL12, 0.628.
(B) Calibration of the five-gene logistic model, with apparent, bootstrap bias-corrected, and ideal curves. (C) Exploratory decision curve. (D) Logistic-regression nomogram based on the five RF-selected genes.
The analysis uses 110 samples (55 per analysis group); calibration used 1,000 bootstrap resamples. Gene selection was not repeated within the calibration bootstrap. These internal summaries do not establish prospective clinical benefit or external validity.Four plots of a five gene expression model with ROC curves, calibration, decision curve and nomogram.A) Line graph of single gene receiver operating characteristic curves.
The y axis label is Sensitivity, unit not shown, ranging from 0.0 to 1.0. The x axis label is 1 minus Specificity, unit not shown, ranging from 0.0 to 1.0 with tick labels 0.0, 0.2, 0.4, 0.6, 0.8, 1.0. Five step like curves are plotted along with a diagonal reference line from 0.0: 0.0 to 1.0: 1.0.
Legend text lists: CCR1, AUC equals 0.615; TRMT2A, AUC equals 0.609; FADS2, AUC equals 0.623; FBXL14, AUC equals 0.651; TTLL12, AUC equals 0.628. The curves rise from near 0.0: 0.0 toward 1.0: 1.0 with multiple small horizontal and vertical steps. Around 1 minus Specificity equals 0.2, the curves are near Sensitivity about 0.2 to 0.3.
Around 1 minus Specificity equals 0.4, the curves are near Sensitivity about 0.55 to 0.65. Around 1 minus Specificity equals 0.6, the curves are near Sensitivity about 0.8 to 0.9. Near 1 minus Specificity equals 1.0, all curves approach Sensitivity 1.0.
B) Line graph of calibration. The x axis label is Predicted probability, unit not shown, with tick labels 0.0, 0.2, 0.4, 0.6, 0.8. The y axis label is Actual probability, unit not shown, with tick labels 0.0, 0.2, 0.4, 0.6, 0.8.
Three lines are shown and a rug of small tick marks appears along the top border. Legend entries read Apparent, Bias corrected, Ideal. The Ideal line is a diagonal from near 0.0: 0.0 toward near 0.9: 0.9.
The Apparent and Bias corrected curves start near the origin, rise above the diagonal around predicted probability near 0.2 to 0.5, then approach the diagonal around predicted probability near 0.6 and end below the diagonal at higher predicted probability near 0.8 to 0.9, where actual probability is around 0.75 to 0.8. C) Line graph of decision curve analysis. The x axis label is Threshold probability, unit not shown, ranging from 0.0 to 1.0 with tick labels 0.0, 0.2, 0.4, 0.6, 0.8, 1.0.
The y axis label is Net Benefit, unit not shown, ranging from 0.0 to 0.5 with tick labels 0.0, 0.1, 0.2, 0.3, 0.4, 0.5. A secondary bottom scale shows cost benefit ratio labels 1:100, 1:4, 2:3, 3:2, 4:1, 100:1. Legend entries read Model, All, None.
The Model curve begins near net benefit 0.5 at threshold probability 0.0, declines gradually to about 0.3 by threshold probability around 0.3, continues downward to around 0.2 by about 0.5, then drops toward near 0.0 around threshold probability about 0.7, dips slightly below 0.0 near about 0.75 and then fluctuates slightly above 0.0 from about 0.8 to 1.0. The All curve slopes downward from near 0.5 at threshold probability 0.0 and crosses near 0.0 around threshold probability about 0.6. The None line lies at net benefit 0.0 across the full threshold probability range.
D) Nomogram style multi scale plot for logistic regression. A top scale is labeled Points with tick labels from 0 to 100 in increments of 10. Below are five gene rows labeled CCR1, TRMT2A, FADS2, FBXL14, TTLL12, each with its own horizontal scale and numeric tick labels.
CCR1 scale shows 10.2, 9.8, 9.4, 9.0, 8.8, 8.4, 8.0, 7.8. TRMT2A scale shows 6.9, 7.1, 7.3, 7.5, 7.7, 7.9, 8.1, 8.3. FADS2 scale shows 7.2, 7.0, 6.8, 6.6, 6.4, 6.2, 6.0, 5.8, 5.6, 5.4.
FBXL14 scale shows 8.7, 8.5, 8.3, 8.1, 7.9, 7.7, 7.5, 7.3. TTLL12 scale shows 7.3, 7.5, 7.7, 7.9, 8.1, 8.3, 8.5. Beneath these is a scale labeled Total Points with tick labels 0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240.
The bottom scale is labeled Risk of Disease with tick labels 0.1, 0.3, 0.7, 0.9, unit not shown.Figure 8Conceptual illustration of the study findings and the proposed biological context. The anatomical elements depict the intestinal barrier, circulating signals, immune cells, and bone-remodeling compartments as an illustrative framework; they do not represent additional measurements. Multi-omics prioritization identified CCR1, TRMT2A, FADS2, FBXL14, and TTLL12, and RANKL-treated RAW264.7 cells showed TRAP-positive osteoclastogenesis accompanied by increased expression of all five proteins.
Solid arrows denote study-supported analytical or experimental relationships, whereas dashed arrows denote proposed relationships. The gut-immune-metabolite-bone connection is a working hypothesis that requires independent mechanistic and participant-level validation.Conceptual biomedical diagram linking gut microbiota, immune cells and bone remodeling with five genes.
The fixed-five-gene ridge assessment included seven datasets and 165 nonoverlapping GSM records. For the primary specification, external AUCs ranged from 0.280 to 0.767. GSE56116 (10 osteoporosis cases and three controls) yielded an AUC of 0.767 (95% CI 0.400–1.000; P = 0.217), whereas GSE56814 yielded 0.392 (95% CI 0.262–0.529).
Across all six specifications, GSE7158 AUCs were 0.601–0.685, providing a limited directional signal in a peak-bone-mass population, but all confidence intervals included 0.5. None of the 42 comparisons reached nominal P < 0.05 (minimum Benjamini-Hochberg q = 0.693). Primary-model Brier scores were worse than the within-cohort prevalence reference in every dataset.
Complete cohort-level results are provided in Supplementary Tables S21–S23 and Figure S5.
Protein Expression During RANKL-Induced Osteoclastogenesis
RANKL treatment induced TRAP-positive multinucleated cells in RAW264.7 cultures (Figure 9A and B). The counts were 0.67 ± 0.58 in controls and 316.67 ± 12.66 in the RANKL-treated group (mean ± SD; n = 3 independent experiments per group; P = 0.000001720). Protein expression was higher in the RANKL-treated group for CCR1 (P = 0.01358), FADS2 (P = 0.04348), FBXL14 (P = 0.001665), TRMT2A (P = 0.02931), and TTLL12 (P = 0.02578) (Figure 9C–H).
These unadjusted, two-tailed Student’s t-test P values were recalculated from the experiment-level values for each outcome. These increases link the five candidate proteins to the differentiation response. Post hoc Benjamini-Hochberg adjustment across the five protein outcomes retained all five associations under the original Student t tests (q = 0.0083–0.0435).
In Welch-test sensitivity analyses, CCR1, TRMT2A, FBXL14 and TTLL12 remained below q = 0.05, whereas FADS2 did not (P = q = 0.0587; Supplementary Table S19). Experiment-level observations and statistical summaries are provided in Supplementary Table S2.Figure 9Protein expression and TRAP-positive cells during RANKL-induced osteoclastogenesis. (A) Representative TRAP staining of control and RANKL-treated RAW264.7 cells; scale bar, 100 μm.
(B) TRAP-positive multinucleated-cell counts (P = 0.000001720). (C) Western blot images with GAPDH as loading control. (D–H) Relative protein expression of CCR1 (P = 0.01358), FADS2 (P = 0.04348), FBXL14 (P = 0.001665), TRMT2A (P = 0.02931), and TTLL12 (P = 0.02578), respectively.
Quantitative data are summarized as mean ± SD with n = 3 independent experiments per group. Each plotted TRAP value represents one independent experiment. P values were recalculated from the experiment-level values using unpaired, two-tailed Student’s t-tests (df = 4) and are unadjusted.
The artwork annotations correspond to *P < 0.05, **P < 0.01, and ****P < 0.0001. Panels A and C show representative images; inference is based on the experiment-level quantitative data. Three parallel wells per group were used in each experiment, with five randomly selected, non-overlapping fields per well assessed for TRAP staining.
Wells and fields were subsamples and were not counted as independent biological replicates. Each experiment-level TRAP value was obtained by averaging counts across the sampled fields and parallel wells.A composite scientific figure of TRAP staining, western blot and bar charts for Control versus RANKL plus.A) Two TRAP-stained cell micrographs labeled, Control and RANKL plus. Control shows sparse stained cells.
RANKL plus shows many stained cells including larger multinucleated cells. A scale bar label reads, 100 micrometer. B) Bar chart of osteoclast counts.
X-axis labels: Control, RANKL plus. Y-axis label: Cell numbers of Osteoclast. Y-axis ticks at 0, 100, 200, 300, 400.
Control shows three plotted points near 0. RANKL plus bar is near 320 with three plotted points around 310 to 330. A significance bracket above the groups is labeled with four asterisks.
C) Western blot with two group headers, Control and RANKL plus, each with three lanes. Rows labeled: CCR1, FADS2, FBXL14, TRMT2A, TTLL12, GAPDH. Molecular weight labels: 41 kilodalton for CCR1, 52 kilodalton for FADS2, 46 kilodalton for FBXL14, 69 kilodalton for TRMT2A, 74 kilodalton for TTLL12, 37 kilodalton for GAPDH.
D) Bar chart of CCR1. X-axis labels: Control, RANKL plus. Y-axis label: Relative expression of CCR1.
Y-axis ticks at 0, 1, 2, 3, 4. Control bar is near 1 with three plotted points around 1.0 to 1.3. RANKL plus bar is near 2.2 with three plotted points around 2.0 to 2.8.
A significance bracket is labeled with one asterisk. E) Bar chart of FADS2. X-axis labels: Control, RANKL plus.
Y-axis label: Relative expression of FADS2. Y-axis ticks at 0, 1, 2, 3. Control bar is near 1 with three plotted points around 1.0 to 1.3.
RANKL plus bar is near 2.0 with three plotted points around 1.7 to 2.6. A significance bracket is labeled with one asterisk. F) Bar chart of FBXL14.
X-axis labels: Control, RANKL plus. Y-axis label: Relative expression of FBXL14. Y-axis ticks at 0, 1, 2, 3, 4.
Control bar is near 1 with three plotted points around 0.9 to 1.2. RANKL plus bar is near 2.8 with three plotted points around 2.6 to 3.3. A significance bracket is labeled with two asterisks.
G) Bar chart of TRMT2A. X-axis labels: Control, RANKL plus. Y-axis label: Relative expression of TRMT2A.
Y-axis ticks at 0, 1, 2, 3, 4. Control bar is near 1 with three plotted points around 1.0 to 1.4. RANKL plus bar is near 2.3 with three plotted points around 2.0 to 3.0.
A significance bracket is labeled with one asterisk. H) Bar chart of TTLL12. X-axis labels: Control, RANKL plus.
Y-axis label: Relative expression of TTLL12. Y-axis ticks at 0, 0.5, 1.0, 1.5, 2.0, 2.5. Control bar is near 1.0 with three plotted points around 1.0 to 1.5.
RANKL plus bar is near 2.0 with three plotted points around 2.0 to 2.4. A significance bracket is labeled with one asterisk.
Discussion
This study combines genetic association analyses, transcriptomic candidate prioritization, exploratory classification, and protein measurements during osteoclast differentiation. The integrated approach identified CCR1, TRMT2A, FADS2, FBXL14, and TTLL12 as candidates and connected their computational prioritization with increased protein expression in RANKL-treated RAW264.7 cells. The study also retained a metabolite association after multiple-testing correction.
Together, these findings provide a focused set of molecular candidates and a reproducible basis for subsequent osteoporosis research.
Microbiota-dependent effects on bone are biologically plausible.7–9 Wang et al reported that Prevotella histicola treatment was associated with preservation of bone in ovariectomized mice and modulation of gut permeability.43 Together with experimental findings for short-chain fatty acids and butyrate, these observations support studying intestinal and immune mechanisms in bone remodeling.4,5 The five candidates offer defined gene-expression and functional readouts for future experiments examining these microbiota-related mechanisms.
CCR1 biology involves both bone-forming and bone-resorbing cells. Hoshino et al reported osteopenia in Ccr1-deficient mice with impaired osteoclast and osteoblast functions.44 This argues against assuming that reduced CCR1 activity must protect bone. CCL2 (MCP-1) signals principally through CCR2; the CCL2-CCR2 pathway should be distinguished from CCR1.
Binder et al linked CCR2-dependent signaling to osteoclast behavior and estrogen-deficiency-associated bone loss.45 These observations place the increased CCR1 expression in the context of receptor-specific contributions to bone remodeling.
Potential osteoblast effects require separate consideration. In addition to the CCR1 findings, descriptive periosteal-cell transplantation work has examined FBXL14 during bone regeneration.46 The present experiments focused on osteoclast differentiation and did not evaluate osteoblast function. Parallel studies of differentiation and mineralization in osteoblasts, including FADS2, TRMT2A, and TTLL12, would help determine how the candidates contribute to both sides of bone remodeling.
The supplementary analyses evaluated the fixed five genes across expression representations, internal resampling procedures, and external cohorts. Regularized logistic regression showed moderate internal discrimination, with mean AUCs of 0.756 for within-sample percentiles and 0.753 for training-reference quantile normalization. These results provide a defined basis for further evaluation of the same five candidates.
The original RF performance estimate describes internal discrimination in the nine-gene analysis. Preprocessing and candidate screening preceded the split, and test-partition importance informed five-gene selection; these steps may make the original performance estimate optimistic.39,40 The separate nomogram and fixed-panel sensitivity analyses remain conditional on that selection. Training-fold reanalysis and external assessment did not establish generalizable prediction.
Differences in cell source, phenotype definition, age, and menopausal status also affect cross-cohort interpretation, and joint ComBat cannot remove these biological differences. Accordingly, LM22 outputs are interpreted as expression signatures rather than bone immune infiltration. The deposited upstream preprocessing was retained, and human transcript and murine protein directions were not uniformly concordant (Supplementary Table S20), reflecting the different biological systems assessed.
The expanded external assessment covers seven datasets while preserving the original five-gene panel. GSE56116 and GSE7158 showed cohort-specific discriminatory signals, but estimates were imprecise and the larger GSE56814 cohort did not reproduce them. External discrimination and probability transfer therefore remain unresolved.
Sample size and phenotype differences, together with possible donor overlap between GSE56116 and GSE230665, are documented in the Supplementary Material; these datasets were not pooled as independent replications.
The corrected MR analysis identified a metabolite association after multiplicity adjustment, with concordant sensitivity-estimate directions. Microbial and immune-trait associations remained nominal, as did the SMR candidate screen, and colocalization did not support a shared signal. These analyses contribute to candidate prioritization; the proposed gut-immune-metabolite relationships remain a biological hypothesis.
The cell experiments complement the computational findings by showing increased expression of all five proteins during TRAP-confirmed osteoclast differentiation in three independent experiments per group. Four candidates retained statistical support across both multiplicity-adjusted Student and Welch analyses, while FADS2 was sensitive to the variance assumption. Gene-specific perturbation, exposure-controlled experiments, and independent participant-level evaluation are the next steps for defining mechanism and diagnostic relevance.
Conclusion
Integrating genetic association analyses, public transcriptomic data, and cell experiments prioritized five osteoporosis-associated candidate genes whose protein expression increased during RANKL-induced osteoclastogenesis. The corrected metabolite association and the moderate internal discrimination of the fixed five-gene panel provide complementary evidence for follow-up. The findings define a focused biological framework for future bone-remodeling studies; predictive use of the panel requires independent participant-level validation.