Integrative Single-Cell RNA Sequencing and Machine Learning Reveals Candidate Plasma Protein-Associated Gene Signatures for Osteoporosis: A Preliminary Exploratory in Silico Study

Introduction

Osteoporosis (OP) is currently the most common metabolic bone disease, characterized by impaired bone cell metabolism that disrupts the balance between bone resorption and new bone formation, ultimately resulting in skeletal fragility and an increased risk of fracture.¹ It was reported that the global prevalence of OP is currently 19.7%, significantly impacting global health.² Meanwhile, with the increasing aging of the population, OP fractures and their complications will place a significant burden on society. At present, dual-energy X-ray absorptiometry (DXA) is the clinical gold standard for OP diagnosis, and bone mineral density (BMD) is also an important indicator for predicting the risk of OP-related fractures. Furthermore, bone mineral density is also a significant indicator for predicting fracture risk.^3,4 Although studies have reported that factors such as estrogen deficiency, inappropriate glucocorticoid use, diabetes, malnutrition, and excessive alcohol consumption are risk factors for OP and fracture occurrence, most people do not pay sufficient attention to OP. Some people only begin to seek treatment of OP after a fracture occurs.⁵ Therefore, early diagnosis and initial risk stratification of OP are crucial for reducing its incidence and facilitating treatment after onset. Discovering more effective diagnostic and therapeutic targets for OP has become the preferred strategy for its prevention and treatment.

Plasma proteins are considered an ideal source of biomarkers due to the non-invasiveness of their collection and dynamic reflection of physiological and pathological states.⁶ Studies have shown that plasma proteins play an important predictive role in the development and progression of various chronic diseases and aging, including Alzheimer’s disease (AD),⁷ cardiovascular and cerebrovascular diseases⁸ and OP. In recent years, multiple studies have highlighted the substantial impact of plasma proteins on bone loss.^9,10 Mendelian randomization (MR) studies have revealed a causal association between 67 plasma proteins and bone mineral density (BMD), with the majority of them influencing OP risk by participating in the formation of the extracellular matrix (ECM).¹¹ Serum proteomics studies have also indicated that 53 proteins, including apolipoproteins and complement components, are associated with OP, suggesting that the serum proteome may serve as a key indicator for assessing skeletal aging.⁹ However, it is currently unclear which plasma protein factors can serve as reliable indicators for OP. Therefore, exploring plasma protein factors with potential regulatory roles holds promise for providing new insights into the diagnosis and treatment of OP. Crucially, we must differentiate between experimentally measured plasma proteins (reflecting direct circulating abundance) and computationally inferred plasma protein-related genes (PRGs). Identified via GWAS, Mendelian randomization, or transcriptomics, PRGs act as genetic or transcriptional proxies rather than direct substitutes for protein levels. Recognizing this distinction is vital for correctly interpreting our integrated findings.^12,13

Various factors, including estrogen deficiency and oxidative stress, act on different cell types such as osteoblasts, osteoclasts, and osteocytes, triggering complex cell-cell interactions that lead to the development of OP.^14,15 Therefore, clarifying the roles of different cell types in the pathogenesis of OP is crucial. Beyond genetic susceptibility and pharmacological interventions, physical activity and rehabilitation training are established modulators of bone metabolism and musculoskeletal integrity.¹⁶ Exercise-induced mechanical loading can influence circulating bone turnover markers, cytokine profiles, and osteogenic signaling pathways, thereby promoting skeletal adaptation and reducing fracture risk.¹⁷ For instance, different exercise modalities have been shown to exert differential effects on pain relief, muscle strength, and functional recovery in patients with degenerative joint diseases such as knee osteoarthritis,¹⁸ underscoring the broader biological link between mechanical stimulation and tissue remodeling. Despite the recognition of these effects, the specific molecular networks through which rehabilitation-related signals regulate bone metabolism at the transcriptomic and single-cell levels remain largely unexplored. This knowledge gap further justifies our integrative approach—namely, the identification of robust and biologically meaningful biomarkers that may also inform future adjunctive therapeutic strategies for osteoporosis based on exercise or rehabilitation.

Single-cell RNA sequencing (scRNA-seq) can precisely profile the gene expression of each individual cell, thereby revealing the complex interactions between different cell types and their unique gene expression signatures.¹⁹ Consequently, integrating bulk RNA-seq, scRNA-seq, and plasma protein data analysis can better combine transcriptomic analysis at the tissue and cellular levels to identify key mechanisms in the onset and progression of OP.

Therefore, this study first successfully identified plasma protein-related biomarkers in OP using bioinformatics methods. Furthermore, enrichment analysis, immune infiltration analysis, and drug prediction were performed to elucidate the potential mechanisms of these biomarkers in OP pathogenesis. Simultaneously, single-cell data analysis was used to characterize the expression patterns of the aforementioned biomarkers in key cell types, providing a basis for in-depth exploration of the functional characteristics of individual cells and revealing heterogeneity within cell populations. Finally, mRNA and Elisa analysis were used to verify whether the expression of biomarkers was consistent with the bioinformatics results. Through systematic exploration of plasma protein genes and their potential mechanisms in OP, this study aims to provide new theoretical support for clinical diagnosis and optimization of treatment strategies for OP.

Materials and Methods

Data Selection

The transcriptome data from RNA-seq datasets (GSE56815 and GSE2208) and the scRNA-seq dataset GSE147287, all of which were linked to OP, were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Among them, GSE56815 (platform: GPL96), designated as the training set, comprised 40 low BMD samples (serving as the OP group) and 40 high BMD samples (serving as the control group). For GSE2208 (platform: GPL96), the validation set, gene expression data were downloaded from 9 low BMD samples versus 10 high BMD samples. Meanwhile, the scRNA-seq dataset GSE147287 (platform: GPL24676) contains one sample of bone marrow obtained from a patient with OP. Additionally, in the Proteome-Phenome-Atlas database (https://proteome-phenome-atlas.com/), Cox regression analysis was performed with sensitivity statistics to mitigate the confounding effects of comorbidities, age, gender, and other covariates. A total of 843 plasma protein-related genes (PRGs) (Table S1) were identified in the general population with OP (p < 0.05).

scRNA-Seq Data Processing

In all samples of GSE147287, to comprehensively investigate the cellular distribution and potential functions across distinct cell populations in OP, an exploratory analysis of scRNA-seq data was performed using the Seurat package (v 5.1.0).²⁰ Raw sequencing data were imported and processed into Seurat objects for downstream analysis. Cells were subjected to stringent quality control with the following inclusion criteria: (1) cells expressing 200–5000 genes; (2) cells with 500–30,000 unique molecular identifiers (UMIs); (3) cells with a mitochondrial gene expression ratio of < 20%. Additionally, genes expressed in fewer than 3 cells were filtered out to minimize noise. After quality control, gene expression values were normalized. And the FindVariableFeatures function was then applied to identify the top 2000 highly variable genes (HVGs). Subsequently, principal component analysis (PCA) was conducted to reduce the dimensionality of the scRNA-seq data. Statistically significant principal components (PCs) were identified using the JackStraw, JackStrawPlot, and ElbowPlot functions (p < 0.05). Clustering analysis was performed using the FindNeighbors and FindClusters functions (resolution = 0.2). To visualize the global clustering landscape of the samples, t-SNE was implemented on the selected PCs. Cell type annotation was performed by referencing canonical marker genes reported in previous studies.^21–23 Osteoblasts differentiated from bone marrow mesenchymal stem cells (BM-MSCs), and the impairment of this differentiation process constitutes a core pathological event in OP.^24,25 To further dissect the cellular heterogeneity of BM-MSCs (as key cells), dimensionality reduction, re-clustering (resolution = 0.1), and re-annotation²¹ were performed. Osteoblasts, a critical cell type derived from BM-MSCs, were the focus of analysis, and genes associated with this cell type were identified.

High-Dimensional Weighted Gene Co-Expression Network Analysis (hdWGCNA)

Based on single-cell level data, the hdWGCNA package (v 0.4.5)²⁶ was utilized, and the MetacellsByGroups function was employed to generate a metacell gene expression matrix for constructing metacells specific to the osteoblasts of the samples. Gene expression profiles associated with the osteoblast-specific module were explored through eigengene connectivity (kME). The ModuleExprScore function of the AUCell algorithm (v 1.22.0)²⁷ was applied to calculate the gene signature score for each module, with genes meeting the criteria of |correlation coefficient (cor) with the module eigengene| > 0.3 and p < 0.05 retained. Genes within the modules that exhibited a strong correlation with osteoblasts were defined as osteoblast-related module genes.

Identification of Candidate Genes

The limma package (v 3.58.1)²⁸ was employed to identify differentially expressed genes (DEGs) between low BMD samples and high BMD samples in GSE56815 (|log₂ fold change (FC)| > 0.1 and p < 0.05).²⁹ The ggvenn package (v 0.1.10)³⁰ was used to identify the overlapping genes among PRGs, osteoblast-related module genes, and DEGs, which were defined as differentially expressed PRGs (DE-PRGs). To elucidate the biological functions of DE-PRGs in OP, clusterProfiler package (v 4.8.1)³¹ was used to perform GO and KEGG analysis (p < 0.05). Furthermore, the DE-PRGs were imported into the STRING database (http://string-db.org, confidence ≥ 0.15, species = Homo sapiens). The protein-protein interaction (PPI) network was visualized using the circlize package (v 0.4.15).³² Subsequently, the MCC, MNC, and Degree algorithms embedded in the Cytohubba plugin of Cytoscape software (v 3.8.2)³³ were employed to identify the top 10 genes, respectively. The overlapping genes among the three sets were defined as candidate genes.

Machine Learning Algorithms

Based on the candidate genes derived from GSE56815, three machine learning algorithms were employed to identify feature genes associated with plasma proteins. LASSO regression was employed to select features (seed = 1), with the optimal regularization parameter identified via 5-fold cross-validation. Then, the glmnet package (v4.1.8)³⁴ was used to train the model with binary logistic regression and L1 regularization. For SVM-RFE (seed = 131, 5-fold cross-validation), the least important features were removed per iteration, and the number of removed features was halved in each subsequent iteration. Error rate curves were generated for all feature subsets, and the optimal subset was selected using the caret package (v 6.0–94).³⁵ A RF model was constructed using the randomForest package (v4.7–1.1)³⁶ (seed = 387). A 5-fold cross-validation procedure was carried out, with the number of trees (ntree) varying from 1 to 100. Feature genes were selected based on the tree count that yielded the minimum error rate. By integrating the results of these algorithms using the ggvenn package (v 0.1.10), the final feature genes associated with plasma proteins were identified.

Identification of Biomarkers and Creation of Nomogram

The expression of the final feature genes was evaluated in GSE56815 and GSE2208 using the Wilcoxon test (p < 0.05), and the feature genes that showed significant differences and consistent expression trends between low BMD samples and high BMD samples were defined as biomarkers. To further appraise the overall predictive capability of biomarkers for OP, a nomogram was established using the rms package (v 6.8–1)³⁷ in GSE56815. Subsequently, calibration curves were generated to exhibit the diagnostic efficacy via rms package (v 6.8–1). The Hosmer-Lemeshow (HL) test served as an indicator for model fitting (p > 0.05). The rmda package (v 1.6)³⁸ was employed to perform decision curve analysis (DCA) to value the clinical utility. The pROC package (v 1.18.5)³⁹ was applied to validate the nomogram’s accuracy. An area under the curve (AUC) value greater than 0.7 was defined as excellent performance.

Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA)

To further investigate the functions and regulatory mechanisms of the biomarkers, GSEA was performed using the whole-genome expression matrix of the GSE56815 dataset. The Spearman correlation analysis was conducted⁴⁰ to calculate the cors between the biomarkers and all other genes via the psych package (v 2.4.3), respectively. Based on the ranking of these cors, GSEA was conducted using the clusterProfiler package (v 4.2.2) (p < 0.05, FDR < 0.25, and |NES)| > 1).

To assess the enrichment status between low BMD samples and high BMD samples, the GSVA package (v 1.50.0)⁴¹ was applied to score KEGG pathways. In parallel, the limma package (v 3.58.1) was utilized to compare the differences in GSVA scores across the two sample groups (|t| > 2 and p < 0.05). Notably, the background gene set (c2.cp.kegg_legacy.v2025.1.Hs.entrez.gmt) was obtained from the MSigDB (http://software.broadinstitute.org/gsea/msigdb/index.jsp).

Immune Infiltration

Immune infiltration analysis identifies diverse immune cell populations, predicts therapeutic responses, and facilitates the development of immunotherapeutic strategies by uncovering the complex interactions between immune cells and biomarkers. In the GSE56815 dataset, infiltration scores for 28 immune cell types⁴² were determined using the ssGSEA algorithm from the GSVA package (v 1.50.0). Differential immune cells (DICs) between low BMD samples and high BMD samples were identified based on the infiltration scores (p < 0.05). Subsequently, Spearman correlation analyses were performed to explore three types of relationships (among biomarkers, among DICs, and between biomarkers and DICs) (|cor| > 0.3 and p < 0.05).

Chromosome Localization and GeneMANIA

Gene localization is of great significance for investigating the structure, function, and interactions of biomarkers. The RCircos package (v 1.2.2)⁴³ was applied to visualize the distribution of these biomarkers on chromosomes. The biomarkers were input into the GeneMANIA (http://www.genemania.org) to construct a gene-gene interaction network, which facilitated the identification and annotation of the biological functions of these biomarkers.

Drug Prediction and Molecular Docking

To provide a basis for the clinical application of OP treatment, the DGIdb (https://dgidb.org/) was utilized to explore potential targeted drugs for the identified biomarkers. The top 15 drugs ranked by the Interaction Score (a specific metric of DGIdb) were visualized. Among the predicted drugs targeting the biomarkers, those that had been approved for clinical use and exhibited the highest Interaction Score were selected for molecular docking analysis, aiming to evaluate their binding affinity. The molecular structures of the drugs were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The crystal structures of biomarkers were obtained from the Protein Data Bank (PDB, http://www.rcsb.org/). CB-Dock2 (https://cadd.labshare.cn/cb-dock2/php/index.php) was employed to predict the binding sites and affinities of protein-ligand complexes. Binding energy with a value of ≤ −5 kcal/mol was considered favorable binding affinity. Molecular docking was performed as an in silico computational prediction to assess potential binding affinities. The results are preliminary and do not constitute evidence of biological activity or therapeutic efficacy, which must be validated through biochemical and cellular assays.

Characterization of Key Cells

In GSE147287, the CellChat package (v 1.6.1)⁴⁴ was employed to separately analyze molecular interactions and ligand-receptor pairs among different annotated cells. Subsequently, the monocle package (v 2.28.0)⁴⁵ was applied for pseudo-time analysis to explore the differentiation trajectory of key cells. Variations in the expression levels of biomarkers during the differentiation of key cells were analyzed in GSE147287. Finally, the SCENIC package (v 1.3.1)⁴⁶ was employed, in combination with the GRNBoost2 algorithm and a list of human transcription factors (TFs), to infer TF-target gene regulatory networks. The Regulatory Specificity Score (RSS) of key cellular regulators in GSE147287 was calculated, and TFs with active regulatory roles that could target the biomarkers were selected. The RSS scores of key cellular regulators across different cell subtypes were visualized, with the top 10 TFs with the highest RSS scores labeled.

Experimental Validation by RT-qPCR and Plasma Protein Level

The clinic data, lymphocytes and serum of OP patients and controls were collected at Zhejiang Provincial People’s Hospital and informed consent was obtained from all people prior to study commencement. This study was approved by the Ethics Committee of Zhejiang Provincial People’s Hospital [approval number ZJPPHEC 2026O (057)]. Venous blood (ethylenediaminetetraacetic acid-anticoagulated) was drawn from patients with OP/osteopenia (n=6) and control (n=6), each group consisted of three males and three females, and lymphocytes were isolated. Total RNA of lymphocytes was extracted by using Trizol method (Invitrogen, United States). The total mRNA (1ug) was reverse transcribed into cDNA (Takara, Japan). And quantitative real time PCR was conducted as previously described⁴⁷ by LightCycler480 (Roche, United States). The primers were as used as previously described.^48,49 The primers were as follows: CTSD: forward primer 5′-GCAAACT GCTGGACATCGCTTG-3′, reverse primer 5′-GCCATAGTGGATGTCAAACGAGG-3′, glyceraldehyde 3-phosphate dehydrogenase (GAPDH): forward primer 5′-CATCATCCCTGCCT CTACTG-3′, reverse primer 5′-GCCTGCTTCACCACCTTC-3′, JUN: forward primer 5′-CTTTTTCGGCACTTGGAGG-3′, reverse primer 5′-GTCCGAGAGCGGACCTTATG-3′ (Sangon Biotech, Shanghai, China).

Further, venous blood was collected from individuals aged from 60 to 65 with normal bone mineral density (n=10) and those with osteopenia/OP (n=10) into pre-chilled tubes and centrifuged (1000 g for 10 min at 4°C) and serum was frozen in liquid nitrogen and stored at −80°C until analyses. And the clinical characteristics of patients was showed in Table S10. An enzyme-linked immunosorbent assay (ELISA) kit for cathepsin D (Wuhan USCN, Wuhan, China) was used for the assay according to the manufacturer’s instructions. Briefly, plasma samples were diluted with Dulbecco’s phosphate-buffered saline (DPBS), and samples and standards were incubated with antibody-coated wells for 1h. After discarding the solutions, the wells were incubated with detection reagent A for 1 h, followed with three times of wash. Then the plate was incubated with detection reagent B for 30mins, followed with three times of wash. Finally, the reaction was then stopped by adding stop solution, and the plate was read immediately at a wavelength of 450 nm using a plate reader (MQX200, Bio-Tek, USA).

Statistical Analysis

All data were processed using R software (v 4.3.3) and GraphPad Prism (v 10). In addition, the Wilcoxon test and the t-test were employed in this study to assess intergroup differences, with the significance threshold set at p < 0.05. The overall analytical workflow for this study is shown in Figure 1.

Figure 1 Analysis flowchart.

Results

Cell Clusters Identification in scRNA-Seq Dataset

Initially, after integrating and filtering the original data of GSE147287 for OP, the data contained 49,436 cells and 9654 genes before quality control and 20,435 cells and 7709 genes after quality control (Figure S1A). After normalizing all samples and selecting the top 2000 HVGs (the top 10 HVGs were: GNLY, IBSP, IGKC, IGLC2, FABP4, IGHA1, IGHG1, DEFA3, RNASE1, HMOX1), dimensionality reduction was performed, and the top 20 PCs were selected (p < 0.05, Figure S1B–D). Then, 14 distinct cell clusters were identified and 8 cell types were annotated (Figure 2A and B). Bubble plots were employed to display marker gene expressions in different cell clusters (Figures S1E and 2C). Eight specific cell types were identified, including BM-MSCs (LEPR, NGFR, ENG, THY1, NT5E), neutrophils (FCGR3B, CSF3R, CXCR2, G0S2), monocytes (CTSS, FCN1, LYZ, CD14), T cells (CD3D, CD3G, TRAC), B cells (CD79A, IGHM, IGHA2, CD19), natural killer (NK) cells (NCAM1), nucleated red blood cells (NRBC, HBA1), and plasmacytoid dendritic cells (pDC, IL3RA, CLEC4C). Previous studies have demonstrated that BM-MSCs undergo progressive senescence with advancing age or estrogen deficiency, which is primarily characterized by a differentiation imbalance: enhanced adipogenic differentiation coupled with impaired osteogenic potential, ultimately leading to reduced bone formation.⁵⁰ Additionally, mechanisms including epigenetic alterations, impaired autophagic function, and accumulated oxidative stress have been recognized as the drivers of OP.⁵¹ To further dissect the heterogeneity of BM-MSC subsets, secondary dimensionality reduction and clustering analysis were performed on BM-MSCs (15 PCs, p < 0.05, Figure S2A–C), which identified three distinct cell clusters (Figures 2D and S2D). Three specific cell types were identified, including osteoblasts (COL1A1, ALPL), adipocytes (MGP, APOD), and terminally differentiated cells (NCAM1, OMD) (Figure 2E and F).

Figure 2 Identification of cell clusters and BM-MSC subtypes in osteoporosis scRNA-seq data. (A) t-SNE plot showing 14 cell clusters identified from the integrated scRNA-seq dataset of osteoporosis (resolution = 0.2). (B) Annotation of the 14 clusters into eight major cell types based on canonical marker genes, including BM-MSCs (LEPR, NGFR, ENG, THY1, NT5E), neutrophils (FCGR3B, CSF3R, CXCR2, G0S2), monocytes (CTSS, FCN1, LYZ, CD14), T cells (CD3D, CD3G, TRAC), B cells (CD79A, IGHM, IGHA2, CD19), NK cells (NCAM1), NRBC (HBA1), and pDC (IL3RA, CLEC4C). (C) Bubble plot showing the expression patterns of representative marker genes across the annotated cell types. (D) t-SNE plot showing the re-clustering of BM-MSCs into three subclusters. (E) Annotation of BM‑MSC subclusters into osteoblasts (COL1A1, ALPL), adipocytes (MGP, APOD), and terminally differentiated cells (NCAM1, OMD). (F) Bubble plot showing the expression of representative marker genes in the three BM-MSC-derived subtypes. Abbreviations: BM-MSCs, bone marrow mesenchymal stem cells; scRNA-seq, single-cell RNA sequencing; NK, natural killer; NRBC, nucleated red blood cell; pDC, plasmacytoid dendritic cell.

Identification of Osteoblast-Related Module Genes by hdWGCNA

Osteoblasts differentiate from BM-MSCs, and impaired differentiation capacity of osteoblasts is a core pathological process of OP.^24,25 Meanwhile, osteoblasts can also regulate the progression of OP by modulating multiple pathways such as metabolism, immunity, and oxidative stress. To further explore the key markers of osteoblasts, hdWGCNA analysis was performed with a soft power of 5 (R² = 0.809, average connectivity = 366.458) to construct a gene co-expression network, which successfully identified three distinct gene co-expression modules Figure 3A and B). Module correlation analysis revealed a strong positive correlation between the brown and blue modules Figure 3C), and both modules showed a strong association with osteoblasts (cor = 0.64, p < 0.0001) Figure 3D and E). To investigate the functional roles of genes in these two modules, genes with significant associations with osteoblasts were screened, and 2453 osteoblast-related module genes were obtained (Table S2). These findings not only verified the robustness of the results but also suggested the potential regulatory role of these osteoblast-related module genes in the pathogenesis of OP.

Figure 3 Identification of osteoblast-related gene modules by hdWGCNA. (A) Soft-threshold power analysis for scale-free topology and mean connectivity to determine the optimal soft‑thresholding parameter for constructing a weighted gene co-expression network specifically in osteoblasts. (B) Hierarchical clustering dendrogram of genes, with module assignments shown in color bars below. Three co-expression modules (brown, blue, grey) were identified using dynamic tree-cutting. (C) Correlation matrix among the identified gene modules. (D) Module signature scores (calculated via AUCell algorithm) for each module in osteoblasts. “+” indicates high activity, and “−” indicates low activity. (E) Heatmap showing the correlations between gene modules and BM-MSC-derived subtypes, highlighting the strong associations of the blue and brown modules with osteoblasts. Red indicates a positive correlation, and blue indicates a negative correlation. The darker the color, the stronger the correlation. Abbreviations: hdWGCNA, high-dimensional weighted gene co-expression network analysis; BM-MSCs, bone marrow mesenchymal stem cells.

Identification and Functional Evaluation of DE-PRGs

In GSE56815, 1214 DEGs (Table S3) were identified between distinct BMD samples. Among them, 648 genes were upregulated and 566 genes were downregulated in the low BMD samples (Figure 4A and B). The intersection of 843 PRGs, 2453 osteoblast-related module genes, and 1214 DEGs was used to identify 15 DE-PRGs, namely EGFR, CTSD, EIF2S2, PTX3, CD99, JUN, ENPP2, ATRAID, HLA-E, BSG, ADAM15, CANT1, ACTA2, LGALS3BP, and CTSO (Figure 4C). A total of 343 GO terms enriched for DE-PRGs included symbiotic interaction, entry into host cell, vesicle lumen, tertiary granule lumen, polysaccharide binding, and exogenous protein binding (p < 0.05, Table S4). KEGG demonstrated that DE-PRGs were enriched in 33 KEGG pathways, such as relaxin signaling pathway, apoptosis pathway, and estrogen signaling pathway (p < 0.05, Figure 4D and Table S5). Based on the STRING database, the protein interactions among 14 of the 15 DE-PRGs (one gene lacked interaction data) involving 31 paired interactions were investigated. Within this network, EGFR, JUN, CTSD, BSG, and LGALS3BP exhibited extensive interactions with other proteins, thereby playing a core regulatory role in the biological network (Figure 4E). The networks of MCC, MNC, and Degree algorithms all consisted of 10 nodes and 25 edges. The top 10 important candidate genes identified from these networks were consistently CD99, JUN, LGALS3BP, PTX3, BSG, CTSD, ACTA2, HLA-E, EGFR, and ENPP2 (Figure 4F).

Figure 4 Identification and functional characterization of differentially expressed plasma protein-related genes in osteoporosis. (A) Volcano plot showing differentially expressed genes between low-BMD and high-BMD samples in GSE56815. Each dot in the plot represents a single gene. Red dots indicate genes with upregulated expression, while blue dots denote genes with downregulated expression. The horizontal axis corresponds to log2FC, the log-transformed value of fold change in gene expression level. The vertical axis is defined as −log10(p value), the negative base-10 logarithm of the P value. (B) Circular heatmap showing the expression patterns of differentially expressed plasma protein-related genes (DE-PRGs). (C) Venn diagram showing the overlap among 843 PRGs, 2453 osteoblast-related module genes, and 1214 DEGs, yielding 15 DE-PRGs. (D) GO (biological process, cellular component, molecular function) and KEGG pathway enrichment results for DE-PRGs (p < 0.05), highlighting pathways such as relaxin signaling, apoptosis, and estrogen signaling. (E) PPI network of 14 DE-PRGs (one gene lacked interaction data) constructed via STRING (confidence ≥ 0.15), with edge width representing combined score. (F) Hub gene network from the intersection of the top 10 genes ranked by MCC, MNC, and Degree algorithms in CytoHubba, consistently identifying CD99, JUN, LGALS3BP, PTX3, BSG, CTSD, ACTA2, HLA-E, EGFR, and ENPP2 as candidate genes. Abbreviations: BMD, bone mineral density; DE-PRGs, differentially expressed plasma protein-related genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene; PPI, protein‑protein interaction.

Identification of JUN and CTSD as Biomarkers for OP

In the LASSO algorithm, the model was constructed using the standard error λ value as the minimum criterion (log (lambda.min) = −3.462) to avoid overfitting and improve prediction accuracy, leading to the identification of 7 LASSO-feature genes, namely EGFR, JUN, CTSD, HLA-E, PTX3, CD99, and ENPP2 (Figure 5A). In the SVM-RFE analysis, when the top 7 features were sufficient to support the optimal discriminative performance of the model, the prediction precision of the model reached its highest value for the first time. This process ultimately identified 7 SVM-RFE-feature genes, including EGFR, PTX3, CTSD, LGALS3BP, ENPP2, HLA-E, and JUN (Figure 5B). The minimum error rate was used as the criterion to select the optimal tree threshold of 19. The constructed RF model identified the top 5 RF-feature genes ranked by importance score, namely EGFR, CD99, CTSD, LGALS3BP, and JUN (Figure 5C). The feature genes from the three machine learning methods were intersected to obtain three final feature genes, including EGFR, CTSD, and JUN (Figure 5D). JUN and CTSD exhibited significantly upregulated expression in the high BMD samples (p < 0.05, Figure 5E and F). Collectively, these results identified JUN and CTSD as biomarkers associated with OP.

Figure 5 Machine learning-based identification of plasma protein-related biomarkers in osteoporosis. (A) LASSO regression analysis, including the coefficient profiles and deviance plot, for feature gene selection. Top: LASSO coefficient profiles for candidate genes as a function of log(λ). Bottom: 5-fold cross-validation deviance plot; the optimal λ (log λ = −3.462, 1‑se criterion) was selected to minimize overfitting, identifying 7 LASSO-feature genes (EGFR, JUN, CTSD, HLA-E, PTX3, CD99, ENPP2). (B) SVM-RFE analysis showing the relationship between variable number and model performance, with the optimal feature subset identified at seven variables (EGFR, PTX3, CTSD, LGALS3BP, ENPP2, HLA-E, JUN). (C) Left: Random forest error rate as a function of the number of trees (ntree), with the minimum error at ntree = 19. Right: Variable importance plot ranking the top 5 genes (EGFR, CD99, CTSD, LGALS3BP, JUN) by importance score. (D) Venn diagram showing the overlap of feature genes (EGFR, CTSD, and JUN) identified by LASSO, SVM-RFE, and random forest. (E and F) Violin plots showing the expression of EGFR, CTSD, and JUN in the GSE56815 training set (E) and GSE2208 validation set (F). Only JUN and CTSD showed significantly consistent upregulation in high-BMD samples across both cohorts (Wilcoxon test), thus defined as biomarkers. ns represents insignificant, *p<0.05, **p<0.01, ***p<0.001. Abbreviations: LASSO, least absolute shrinkage and selection operator; SVM-RFE, support vector machine-recursive feature elimination.

An OP diagnostic nomogram model was developed based on the biomarkers (JUN and CTSD) from GSE56815. The higher the total points, the greater the probability of suffering from OP (Figure 6A). The calibration curve nearly overlapped with the ideal prediction curve, indicating that the model had reliable predictive accuracy (p = 0.507) (Figure 6B). The model curve maintained a distinct distance from the two extreme curves, indicative of favorable predictive performance. Moreover, the nomogram yielded a higher net benefit rate compared with any single gene, further corroborating the superior diagnostic efficacy of the established model (Figure 6C). The AUC of the nomogram was 0.762 in the GSE56815, and the AUC of the nomogram was 0.922 in the GSE2208, indicating the reliability of the model’s diagnostic efficacy (Figure 6D and E).

Figure 6 Diagnostic nomogram based on the identified biomarkers. (A) Nomogram integrating JUN and CTSD for individualized prediction of OP probability. Points are assigned to each biomarker based on its expression value; total points correspond to predicted risk, **p<0.01. (B) Calibration curve for the nomogram in GSE56815. The dashed diagonal line represents ideal prediction; the solid curve shows actual performance, with the Hosmer-Lemeshow test p = 0.507 indicating good fit. (C) Decision curve analysis comparing the net clinical benefit of the nomogram (blue) versus individual genes (red and yellow) and the “All” (gray) or “None” (black) strategies. The nomogram provides good net benefit across a wide range of threshold probabilities. (D and E) Receiver operating characteristic curve showing the diagnostic performance of the nomogram in the GSE56815 (D) and GSE2208 (E) datasets.

Biological Mechanisms Associated with Biomarkers in OP

GSEA results demonstrated that JUN was enriched in 10 distinct pathways, with key pathways encompassing the ribosome, proteasome, and spliceosome pathways (p < 0.05, Figure 7A and Table S6). By contrast, CTSD was enriched in 23 pathways, including the lysosome, RNA degradation, and citrate cycle (TCA cycle) pathways (p < 0.05, Figure 7B and Table S7). Notably, both CTSD and JUN were significantly enriched in the five overlapping pathways, such as ribosome, proteasome, lysosome, and spliceosome (p < 0.05, Figure 7C and D). Furthermore, GSVA results revealed that a total of 48 pathways exhibited significant intergroup differences (p < 0.05). The representative pathways significantly upregulated in the low BMD samples included ABC transporters and aminoacyl tRNA biosynthesis (p < 0.05, Figure 7E and Table S8).

Figure 7 Functional pathway analyses of JUN and CTSD in osteoporosis. (A) GSEA enrichment plots for JUN-correlated genes (Spearman correlation with whole-genome expression in GSE56815). Significant pathways include ribosome, proteasome, and spliceosome (p < 0.05, FDR < 0.25, |NES| > 1). (B) GSEA plots showing representative pathways significantly associated with CTSD. (C) Venn diagram showing the overlap of significantly enriched pathways between JUN and CTSD. (D) Heatmap showing the normalized enrichment scores of the overlapping pathways shared by JUN and CTSD, *p<0.05, **p<0.01, ***p<0.001. (E) GSVA plots showing representative KEGG pathways with significant differences between low-BMD and high-BMD samples. *p<0.05, **p<0.01, ***p<0.001. Abbreviations: GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; NES, normalized enrichment score, FDR, false discovery rate.

Immune Landscape in OP

To conduct a more in-depth exploration of the immune microenvironment in OP, immune cell infiltration analysis was carried out using data from the GSE56815 dataset, clarifying the intricate interactions among immune cells. Compared with the low BMD samples, the infiltration levels of activated dendritic cells and mast cells were markedly upregulated in the high BMD samples (p < 0.05, Figure 8A and B). A positive correlation between activated dendritic cells and mast cells was found (cor = 0.61, p < 0.001). All differential immune cells exhibited positive correlations with the biomarkers: JUN showed strong positive correlations with activated dendritic cells (cor = 0.37, p < 0.001) and mast cells (cor = 0.57, p < 0.001), and CTSD also showed positive correlations with these two immune cell types (Figure 8C).

Figure 8 Immune landscape associated with the identified biomarkers in osteoporosis. (A) Heatmap of ssGSEA-based infiltration scores for 28 immune cell types across GSE56815 samples. (B) Boxplots showing activated dendritic cells and mast cells as the two cell types with significantly different infiltration levels between low- and high-BMD groups (Wilcoxon test), ns represents insignificant, *p<0.05. (C) Spearman correlation matrix among JUN, CTSD, activated dendritic cells, and mast cells. JUN showed strong positive correlations with activated dendritic cells (cor = 0.37, p < 0.001) and mast cells (cor = 0.57, p < 0.001); CTSD also correlated positively with both, though weaker. ns represents insignificant, **p<0.01, ***p<0.001. Abbreviation: ssGSEA, single-sample gene set enrichment analysis; BMD, bone mineral density.

Comprehensive Analysis of Biomarkers in OP

Chromosomal localization showed that JUN and CTSD are located on chromosome 1 and chromosome 11, respectively (Figure 9A). To further explore the interactions between the hub genes, the GeneMANIA database was employed. The co-expression network revealed associations between the biomarkers and 20 related genes, such as FOS, CTSA, ATF3, MAPK8, IGFBP7, MAPK9, ATF4, CREB5, IFI30, and KDM6B. These genes are likely involved in the same biological processes, including response to cadmium ion, cellular response to oxidative stress, Fc receptor signaling pathway, and RNA polymerase II transcription regulator complex (Figure 9B). The hub genes were submitted to the DGIdb database, and the results revealed that CTSD was linked to two candidate drugs, among which one (streptozocin) was a clinically approved agent and one (C3TD879) was an investigational compound. In contrast, JUN was associated with 46 drugs, including 24 approved pharmaceuticals (such as bruceantin and ciprofibrate) and 22 unapproved candidates (such as risolidone and sergeolide) (Table S9; only the top 15 drugs ranked by Interaction Score are presented herein) (Figure 9C). Subsequently, the approved drugs with the highest Interaction Score predicted to target these biomarkers were selected for molecular docking validation. The results revealed that streptozocin bound to CTSD with a binding energy of −6.5 kcal/mol, interacting via residues including GLY-79, SER-80, THR-125, ASP-33, and GLY-35 (Figure 9D). Bruceantin bound to JUN with a binding energy of −8.5 kcal/mol, interacting via residues including GLN-669 and ARG-286 (Figure 9E). This indicated that the docking results predicted favorable binding conformations, providing a computational basis for future experimental validation.

Figure 9 Comprehensive characterization of JUN and CTSD and prediction of candidate drugs. (A) Chromosomal localization of JUN and CTSD. (B) GeneMANIA network showing genes functionally associated with JUN and CTSD. (C) Drug-gene interaction network showing the top predicted drugs targeting the biomarkers. (D) Molecular docking model of streptozocin with CTSD. (E) Molecular docking model of bruceantin with JUN.

Characterization of BM-MSCs in OP

The cell-cell communication network was visualized to characterize the intensity and interaction counts of intercellular communication. Notably, the crosstalk between BM-MSCs and NK cells exhibited the highest activity level, ranking first in both interaction counts and communication intensity. These findings suggested that our exploratory analysis raises the possibility that BM-MSCs may influence NK cell activity, although this observation is derived from a single sample and requires confirmation. In contrast, neutrophils displayed minimal interactions with other cell types, implying that they might be in an immunosuppressed or senescent state within the chronic bone metabolic microenvironment (Figure 10A and B). Immune cells could activate BM-MSC functions via pathways such as SPP1, MDK, and FGF, and these pathways might either promote bone repair or exacerbate bone loss. Conversely, BM-MSCs could modulate immune responses and matrix remodeling through key signaling axes including CXCL12-CXCR4 and ANGPTL4-CDH11, thus forming a dynamic feedback loop (Figure 10C). BM-MSCs underwent two critical transition nodes during their development and differentiation. Adipocytes represented an early developmental stage, whereas terminally differentiated cells corresponded to the late stage of BM-MSC development and differentiation. The developmental timeline of macrophages was divided into five stages: Stage 1 was the early developmental stage and Stage 5 was the late developmental stage (Figure 11A). The expression patterns of biomarkers across pseudo-time trajectories were unveiled in BM-MSCs. CTSD was predominantly highly expressed in the early-to-middle phases of differentiation, whereas JUN exhibited robust expression mainly during the middle and terminal phases (Figure 11B). Additionally, this study revealed that BM-MSC subtypes were closely correlated with TFs. Specifically, TFs including EGR1, IRF1, JUN, and JUNB exhibited high activity in the osteoblast subtype, suggesting that these cells were in an intermediate state between osteogenic function execution and pathological response (Figure 11C). Moreover, the top 10 scoring TFs varied across the three distinct key cell subsets, which may indicate the presence of divergent transcriptional regulatory profiles among these cellular subtypes (Figure 11D).

Figure 10 Characterization of BM-MSC communication in osteoporosis. (A) Intercellular communication network (CellChat) showing the number of interactions among all cell types. BM-MSCs and NK cells showed the highest interaction counts; neutrophils had minimal interactions. (B) Heatmaps showing the interaction number and interaction strength among cell types. (C) Bubble plots showing the major incoming and outgoing signaling pathways involved in BM-MSC communication. Abbreviation: BM-MSCs, bone marrow mesenchymal stem cells.

Figure 11 Characterization of BM-MSC differentiation trajectory, and transcriptional regulation in osteoporosis. (A) Monocle2 pseudotime trajectory of BM-MSC differentiation, with cells colored by pseudotime (left), cell state (middle), and subtype (right). Two bifurcation nodes were identified; adipocytes appeared at early stages and terminally differentiated cells at late stages. (B) Dynamic expression patterns of JUN and CTSD along the BM-MSC differentiation trajectory. (C) Heatmap showing the activity of representative transcription factors across osteoblasts, adipocytes, and terminally differentiated cells. “+” indicates that this transcription factor potentially exerts a positive regulatory effect. (D) Regulon specificity score plots showing the top transcription factors in osteoblast, terminally differentiated, and adipocyte subtypes. Abbreviation: BM-MSCs, bone marrow mesenchymal stem cells.

Expression Confirmation of JUN and CTSD

To further study the expression of CTSD and JUN, 4 mL of venous blood (ethylenediaminetetraacetic acid-anticoagulated) was drawn from patients with OP and control, and lymphocytes were isolated, followed by RT-qPCR assays. The results showed that the mRNA expression of CTSD was significantly upregulated compared to the control group, while the expression of JUN showed no significant difference. Furthermore, we selected 20 serum samples from individuals aged from 60 to 65 with normal bone mineral density and those with osteopenia or osteoporosis, the baseline characteristics showed in Table S10. We then used ELISA to detect serum CTSD activity. The results showed that CTSD activity was significantly higher in the osteopenia or osteoporosis group than in the normal bone mineral density group, suggesting that CTSD may be associated with osteoporosis status (Figure 12).

Figure 12 Analysis of the RT-qPCR expression and activity of CTSD and JUN. The relative mRNA expression levels of CTSD and JUN, and the activity of serum CTSD detected by ELISA, **p<0.01.

Discussion

Plasma proteins play a critical role in numerous biological processes and the abnormal alterations of these proteins can lead to changes in physiological conditions. In recent years, there has been increasing research focusing on the relationship between plasma proteins and OP, with multiple studies reporting associations between plasma proteins and BMD as well as OP. Aparicio-Bautista et al⁵² reported that levels of several plasma proteins, including APOA1, SHBG, and FETB, differed in postmenopausal women with OP and osteopenia. Another study also found differential expression of the plasma protein ANXA2 between subjects with low and high BMD,⁵³ suggesting that plasma proteins are important in the development of OP. By integrating whole-transcriptome analysis, single-cell sequencing, and machine learning approaches, this study systematically identified JUN and CTSD as plasma protein-related biomarkers for osteoporosis (OP), which were partially validated by RT-qPCR and ELISA assays. These two biomarkers are enriched in bone metabolism-related pathways and correlate with the infiltration of specific immune cells. Furthermore, the nomogram model constructed based on these biomarkers demonstrated promising diagnostic potential in retrospective datasets. However, all computational findings—particularly those related to immune regulation, molecular docking, and single-cell inference—are exploratory in nature and derived from a limited sample size. Independent validation through functional experiments and larger-scale prospective cohort studies is required before any translational medicine conclusions can be drawn.

Jun proto-oncogene (Jun) is one of the AP-1 transcription factor family proteins and has been reported to be mainly involved in fibrotic diseases and regulate key cellular processes such as the cell cycle.⁵⁴ JUN also plays an important regulatory role in the skeletal system. Lerbs et al demonstrated using different Jun-induced mouse models that JUN can promote osteogenic differentiation and bone formation, as well as aiding fracture repair.⁵⁵ Mechanistically, it is reported that JUN may promote bone formation by regulating osteogenic differentiation of BMSCs through the Hedgehog and Wnt signaling pathway.^55–57 Notably, JUN’s role in bone metabolism is double-edged: Under stimulation by inflammatory signals (such as TNF-α and IL-6), phosphorylation or upregulation of c-Jun can enhance osteoclast activity and promote bone resorption.^58,59 Furthermore, our RT-qPCR validation did not detect significant differences in JUN mRNA expression between OP patients and controls. This observation can be attributed to multiple factors. First, the biological function of JUN is primarily dictated by its protein levels and post-translational modifications (PTMs) rather than mRNA abundance. As a core component of the AP-1 transcription complex, JUN’s transcriptional activity is strictly regulated by JNK-mediated phosphorylation at Ser63 and Ser73, as well as acetylation modifications that affect its stability.^53,54 Therefore, mRNA levels may not directly reflect the functional activity of JUN. Second, JUN expression is subject to post-transcriptional regulation. Studies have shown that the RNA-binding protein ZFP36 suppresses JUN translation by binding to the 3’UTR of JUN mRNA, suggesting that protein expression can be uncoupled from its mRNA abundance.⁶⁰ More broadly, the weak correlation between mRNA and protein levels of transcription factors is a well-recognized phenomenon. Third, JUN exerts a dual regulatory role in bone metabolism: it promotes osteogenic differentiation via the Hedgehog pathway while simultaneously driving osteoclast differentiation through the c-Fos/c-Jun axis.^50–54 Given this functional complexity, measuring JUN mRNA alone is insufficient to comprehensively evaluate its role in OP. Future studies should incorporate protein-level and phosphorylation status analyses to better elucidate the involvement of JUN in this disease.

CTSD is a lysosomal aspartyl protease,^61,62 primarily expressed in humans as an inactive proenzyme form.⁶³ Studies have shown that the CTSD tissue-specific expression plays a crucial role in determining the development of different tissues.⁶⁴ Recent research suggested that CTSD may be involved in the pathogenesis of OP^65,66 and has the potential to be a biomarker for OP.⁶⁵ However, current studies on the relationship between CTSD and OP are limited, and the underlying mechanisms remain unclear.⁶⁷ Animal studies by Sun et al have found that CTSD can influence the osteogenic capacity of rat BMSCs, but the mechanism remained unclear.⁶⁸ Conversely, research by Song et al suggested that CTSD can weaken osteogenic differentiation of BMSCs and new bone formation, which may be related to the ability of CTSD to degrade bone matrix proteins, regulate osteoclast activity, and indirectly affect the expression and function of plasma proteins related to bone metabolism.^67,69 The seemingly contradictory results mentioned above suggest that the role of CTSD in bone metabolism may be context-dependent, and its specific regulatory mechanisms require further investigation. Furthermore, the apparent discrepancy in CTSD mRNA levels between the public transcriptomic data and our experimental validation may stem from multiple factors. First, as a lysosomal aspartic protease, CTSD synthesis requires signal peptide cleavage, glycosylation, and multiple proteolytic cleavages to transition from an inactive proenzyme to a mature enzyme.^56,57 Consequently, mRNA abundance does not always parallel its proteolytic activity. Second, the cellular composition may vary across different datasets; since CTSD expression patterns can differ among immune cell subsets, this could affect the comparability of transcriptomic results. More importantly, the function of CTSD in bone metabolism primarily depends on its enzymatic activity rather than its transcriptional level. CTSD can degrade non-collagenous proteins in the bone matrix and promote osteoclast-mediated bone resorption,⁶³ and animal studies have confirmed that CtsD deficiency leads to a significant decrease in bone mass.^61,62 Therefore, the significantly elevated serum CTSD activity observed in OP patients in our study (via ELISA) provides more direct functional evidence than mRNA trends, supporting the rationale for its use as a candidate biomarker. Notably, regarding the specificity of JUN and CTSD, it is necessary to carefully delineate their value as biomarkers for osteoporosis (OP) from their broader associations with systemic inflammation or aging. As a core member of the AP-1 transcription factor family, JUN promotes inflammatory bone erosion in rheumatoid arthritis (RA) by regulating COX2 expression in macrophages, and the JNK/c-Jun signaling pathway is highly activated in RA synovial fibroblasts.^70,71 However, its function in bone metabolism exhibits a cell type-dependent dual regulatory pattern. In osteoblast precursors, JUN significantly upregulates RUNX2 expression to drive osteogenic differentiation,^55,72 whereas in an inflammatory microenvironment, c-Jun phosphorylation may enhance osteoclast activity.⁷³ This functional duality suggests that the altered expression of JUN in OP must be interpreted within the context of bone-specific cells. Similarly, although CTSD, a lysosomal protease, is found at reduced levels in the plasma of Alzheimer’s disease (AD) patients and has been explored as a potential AD biomarker,⁷⁴ its regulation in bone tissue may be unique. Yan et al demonstrated that CTSD is differentially expressed in the bone tissue of OP patients. Recombinant human CTSD downregulates the expression of osteogenic markers (ALP, RUNX2, and COL1A1), and exogenous CTSD reduces bone mineral density in mice, mimicking the phenotype of the ovariectomized (OVX) model.⁶⁷ Furthermore, CTSD knockout mice also exhibit a significant reduction in bone mass.⁷⁵ Collectively, these findings indicate that despite the involvement of JUN and CTSD in broad physiological and pathological processes, their expression and functions in bone tissue may be tissue-specific. Therefore, their value as OP biomarkers requires validation through bone microenvironment-specific assessments and large-scale cohort studies, rather than being simply attributed to concomitant systemic inflammation or aging.

To further explore the potential mechanisms of CTSD and JUN in OP, we performed pathway enrichment analysis. The results showed that both were significantly enriched in pathways related to ribosomes, proteasomes, lysosomes, spliceosomes, and Parkinson’s disease (PD). Notably, these pathways are closely related to intracellular protein homeostasis, suggesting that CTSD and JUN may show correlation with the pathogenesis of OP by synergistically regulating protein synthesis, folding, degradation, and RNA splicing. Previous studies have reported that various ribosomal proteins can affect bone metabolism, including RPS17, RPS12, RPL31, RPL34, and PRL29.^76–79 These proteins may affect ribosome biogenesis and protein translation speed which can promote the translation of osteoblast differentiation-related factors, while inhibition of ribosome assembly can reduce osteoblast activity.^80,81 The regulatory mechanism involved is currently unclear, but some researchers speculate that it may be related to inflammation activation or the impact of protein translation efficiency relating to bone formation and the bone microenvironment.^76,79 Researchers have proposed that there are some similarities in the pathogenesis of PD and OP.⁸² Although less studied, PD-related pathways have some evidence suggesting that proteins like α-synuclein may influence bone remodeling via the circulatory system or bone marrow microenvironment.⁸³ Furthermore, the improved bone microstructure observed in LRRK2 knockout mice suggests a possible interaction between this pathway and Wnt signaling.⁸⁴

The proteasome and lysosomal pathways are closely related to bone metabolism, and their correlations have been extensively reported in many studies. Currently, several proteasome inhibitors (PIs) have been reported to control the degradation of key proteins in osteoblast differentiation, thereby promoting bone anabolism and stimulating osteoblast differentiation and fracture healing. Furthermore, studies have shown that PIs can also effectively prevent osteocyte death, thus reducing bone tissue loss.⁸⁵ Lysosomes are closely related to various bone cell functions and bone metabolism. Specifically, the process that osteoblasts maintains mineralization rely on lysosome-mediated matrix vesicle transport and ECM component recycling. Osteoclasts inducing bone resorption rely on lysosome-driven acidification and protease release. Osteocytes utilize the lysosome-autophagy axis to maintain the patency of the lacunar-canalicular network and maintain responsiveness to microenvironmental signals.^86–88 Splicing pathways primarily affect bone metabolism by regulating the transcription of key regulatory proteins in osteoblasts or osteoclasts. Existing studies have identified and validated several splice variants associated with key osteoblastogenesis or osteoclastogenesis proteins, such as TAF4, CBFB2, vRANK, and Bcl-xl.⁸⁹ Considering the molecular characteristics of CTSD as a lysosomal protease and JUN as a transcription factor, we hypothesize that CTSD may directly participate in bone matrix degradation and osteoclast activation through the lysosomal pathway, while JUN may regulate the translation and splicing of osteogenesis-related genes through the ribosome/spliceosome pathway. Their co-enrichment in these pathways suggests they may form a functional axis that coordinately regulates the OP process. This finding provides a new direction for further in-depth studies on the specific molecular mechanisms of CTSD and JUN in OP.

The results of immune infiltration analysis in this study suggest a significant difference between the OP group and the control group in the infiltration levels of activated dendritic cells (DCs) and mast cells (MCs). Further analysis indicated that both immune cell types were positively correlated with JUN and CTSD, suggesting that these two biomarkers may participate in the immunopathological process of OP by regulating these cells. DCs are classified as part of the innate immune system. Their primary role is as antigen-presenting cells (APCs), enabling them to recognize pathogens, activate through antigen processing, migrate, present antigens to naive T cells, and secrete cytokines.⁹⁰ In bone metabolism, DCs are closely related to osteoclasts: both originate from hematopoietic stem cells, exhibit highly overlapping gene expression profiles, and share the same surface markers.⁹¹ Under stimulation by inflammatory factors, immature DCs can transdifferentiate into osteoclasts under specific conditions. Conversely, osteoclastic products released from bone destruction areas can also promote the differentiation and activation of immature DCs, forming a feedback loop regulating bone metabolism.⁹² MCs are immune cells within tissues, originating from pluripotent precursor cells in the bone marrow. MCs precursors migrate into tissues where they differentiate and mature.⁹³ Current research suggests that MCs play a bidirectional regulatory role in bone metabolism: on the one hand, activated mature MCs can rapidly release inflammatory mediators such as histamine and IL-6, enhancing osteoclast activity and promoting bone destruction.⁹⁴ On the other hand, MCs can also secrete factors such as TGF-β, IL-12, and IFN-γ, enhancing osteoblast activity and inhibiting osteoclasts, thereby exerting a bone-protective effect.^95,96 This functional duality suggests that the role of MCs in OP may depend on the dynamic balance of local microenvironmental signals. In summary, the positive correlation between JUN and activated DCs and MCs in this study, combined with the critical roles of these two immune cells in bone metabolism, suggests that JUN may participate in the immunomodulatory network of OP by regulating the function of DCs and MCs. In contrast, CTSD exhibits weaker correlations with these immune cells, suggesting its role might be more confined to processes such as lysosome-mediated antigen processing or bone matrix degradation, rather than direct immune regulation. This finding provides a new entry point for a deeper understanding of the molecular mechanisms of JUN in OP.

Through scRNA-seq analysis, this study identifies BM-MSCs as a critical cell type in the pathogenesis of OP. Consistent with other research, the imbalance in osteogenic and adipogenic differentiation of BM-MSCs is considered a core component of OP development.^97–99 This study reveals a dynamic expression pattern of the biomarker JUN during BM-MSC differentiation. We guess that c-Jun overexpression may activates the NF-κB pathway under inflammatory or pathological conditions, inhibiting the activity of the key osteogenic transcription factor RUNX2, thereby impairing osteogenic differentiation and promoting adipogenic differentiation.¹⁰⁰ Targeting and inhibiting c-Jun overactivation can restore the osteogenic capacity of BM-MSCs and improve BMD in OP model mice.¹⁰¹ These findings showed that JUN can be a a potential candidate for therapeutic exploration for OP. However, the scRNA-seq analysis in this study was based on a single OP sample. Consequently, the identified cellular subpopulation dynamics, pseudotime trajectories, and transcription factor regulatory networks are exploratory in nature. The generalizability and mechanistic interpretability of these findings must be validated in larger-scale single-cell cohorts.

In conclusion, as a preliminary exploratory study, this research integrated single-cell RNA sequencing, machine learning, and bulk transcriptomic data to screen and preliminarily validate JUN and CTSD as genetic biomarkers for osteoporosis (OP). Bioinformatic analyses suggested an association between these two genes and bone metabolism-related pathways as well as the immune microenvironment, while single-cell data revealed their expression characteristics during bone marrow mesenchymal stem cell differentiation. However, this study has several important limitations. First, the single-cell sequencing data were derived from only one OP patient. This extremely small sample size fails to represent the cellular heterogeneity of the OP population; thus, the resulting findings, such as cell-cell communication, should be regarded as exploratory. Second, due to the cross-sectional study design, this study can only reveal concurrent associations between variables. It cannot establish causality, nor can it infer the temporal evolution of the disease or predict individual-level risks of disease progression. Third, the current findings remain preliminary. All results are based on computational and correlational evidence, and their clinical significance and therapeutic value require validation through independent cohorts and experimental studies. They should not be overinterpreted as directly translatable clinical diagnostic or therapeutic strategies.

Overall, the primary significance of this study lies in providing novel candidate targets and testable hypotheses for the investigation of the molecular mechanisms underlying OP. Future research urgently requires large-scale, multicenter prospective cohort studies to validate the stability and reproducibility of these biomarkers over a longitudinal timeframe. Concurrently, it is necessary to expand the sample size for single-cell sequencing to include more OP patients and diverse clinical phenotypic subgroups, thereby comprehensively assessing cellular heterogeneity and verifying the generalizability of our findings. Furthermore, in vivo and in vitro functional experiments are essential to elucidate the specific mechanisms of JUN and CTSD in bone metabolism, allowing for a rigorous assessment of their clinical translational potential.

Conclusions

This preliminary exploratory study integrated bulk transcriptomic analysis, single-cell RNA sequencing, and machine learning algorithms to identify JUN and CTSD as plasma protein-related gene biomarkers for OP. Bioinformatic analyses revealed that these two biomarkers are significantly enriched in pathways related to ribosomes, proteasomes, lysosomes, and spliceosomes, suggesting their potential involvement in protein homeostasis regulation during OP pathogenesis. Immune infiltration analysis demonstrated positive correlations between JUN and activated dendritic cells and mast cells. Single-cell characterization further revealed distinct expression dynamics of JUN and CTSD during BM-MSCs differentiation, providing cellular context for their roles in bone metabolism. Experimental validation via RT-qPCR and ELISA confirmed elevated CTSD expression and activity in OP patients. Nevertheless, these findings are derived from computational and correlational evidence with limited sample sizes, and their clinical translational potential requires further validation through large-scale prospective cohorts and functional experiments.