Integrated single-cell and bulk RNA sequencing analysis reveal immune-related biomarkers in postmenopausal osteoporosis

Abstract

Background: Postmenopausal osteoporosis (PMOP) represents as a significant health concern, particularly as the population ages. Currently, there is a paucity of comprehensive descriptions regarding the immunoregulatory mechanisms and early diagnostic biomarkers associated with PMOP. This study aims to examine immune-related differentially expressed genes (IR-DEGs) in the peripheral blood mononuclear cells of PMOP patients to identify immunological patterns and diagnostic biomarkers.

Methods: The GSE56815 dataset from the Gene Expression Omnibus (GEO) database was used as the training group, while the GSE2208 dataset served as the validation group. Initially, differential expression analysis was conducted after data integration to identify IR-DEGs in the peripheral blood mononuclear cells of PMOP. Subsequently, feature selection of these IR-DEGs was performed using RF, SVM-RFE, and LASSO regression models. Additionally, the expression of IR-DEGs in distinct bone marrow cell subtypes was analyzed using single-cell RNA sequencing (scRNA-seq) datasets, allowing the identification of cellular communication patterns within various cell subgroups. Finally, molecular subtypes and diagnostic models for PMOP were constructed based on these selected IR-DEGs. Furthermore, the expression levels of characteristic IR-DEGs were examined in rat osteoporosis (OP) models.

Results: Using machine learning, six IR-DEGs (JUN, HMOX1, CYSLTR2, TNFSF8, IL1R2, and SSTR5) were identified. Subsequently, two molecular subtypes of PMOP (subtype 1 and subtype 2) were established, with subtype 1 exhibiting a higher proportion of M1 macrophage infiltration. Analysis of the scRNA-seq dataset revealed 11 distinct cell clusters. It was noted that JUN was significantly overexpressed in M1 macrophages, while HMOX1 showed a marked elevation in endothelial cells and M2 macrophages. Cell communication results suggested that the PMOP microenvironment features increased interactions among M2 macrophages, CD8⁺ T cells, Tregs, and fibroblasts. The diagnostic model based on these six IR-DEGs demonstrated excellent diagnostic performance (AUC = 0.927). In the OP rat model, the expression of IL1R2 and TNFSF8 were significantly elevated.

Conclusion: JUN, HMOX1, CYSLTR2, TNFSF8, IL1R2, and SSTR5 may serve as promising molecular targets for diagnosing and subtyping patients with PMOP. These results offer novel perspectives on the early diagnosis of PMOP and the advancement of personalized immune-based therapies.

Keywords: Biomarkers; Diagnosis; Immune; Molecular subtype; Postmenopausal osteoporosis.

Fig. 1

The flowchart of this study.

Fig. 2

Identification and Functional Analysis of Immune-Related Differentially Expressed Genes (IR-DEGs). (A) Differential heat map of the top 10 IR-DEGs, with red representing upregulated genes and blue representing downregulated genes. (B) Differential box plot of the top 10 IR-DEGs, with red indicating the PMOP group and blue indicating the control group. (C) The top 10 significant enrichment of KEGG pathways. (D) The top 5 significantly enriched biological functions (BP) revealed by GO functional analysis, the top 4 enriched cellular components, and the top 5 enriched molecular functions. Significance: ∗∗P < 0.01, ∗∗∗P < 0.001.

Fig. 3

Identification of characteristic IR-DEGs using machine learning methods. (A and B) Residuals from the Random Forest (RF) and Support Vector Machine-Recursive Feature Elimination (SVM-RFE) models. (C) Change curves in the predicted true values for each IR-DEG. (D) Change curves in the predicted error values for each IR-DEG. (E) Gene coefficient plots from the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis. (F) Cross-validation curves. (G) Characteristic IR-DEGs obtained by the intersection of genes selected by LASSO and SVM-RFE. (H) Chromosomal positions of the 6 characteristic IR-DEGs.

Fig. 4

Cell clustering analysis and annotation using the tSNE method. (A) Identification of 12 cell clusters based on the tSNE method. (B) Determination of 11 cell types based on cell-specific markers. (C) Specific markers for the 11 cell types. (D) Heatmap illustrating the classification of the 11 cell types.

Fig. 5

Cell cycle and cell communication analysis in PMOP patients. (A) Expression levels of five characteristic IR-DEGs in bone marrow cells. (B) Visualization of differential expression of five characteristic IR-DEGs among different cell types. (C) Proportions of different cell types in various cell cycles. (D) Visualization of the cell cycle distribution among different cell types. (E) Interactions among the 11 different cell types.

Fig. 6

Construction of molecular subtypes in peripheral blood monocytes of PMOP patients based on characteristic IR-DEGs. (A) Heatmap of different K-means clustering matrices (K = 2–4). (B) Consistency clustering scores of different subtypes (K = 2–9). (C) Cumulative distribution function (CDF) plot determining the optimal K value (K = 2). (D) Delta area plot. (E) Principal component analysis (PCA) of PMOP subtypes 1 and 2. (F) Differential expression analysis of the 6 IR-DEGs between Subtypes 1 and 2. Significance: ∗∗P < 0.01, ∗∗∗P < 0.001.

Fig. 7

Analysis of the immune microenvironment in PMOP patients. (A) Heatmap depicting the correlation between the infiltration abundance of 22 immune cell types between the PMOP group and control group. (B) Differential expression analysis of the infiltration abundance of 22 immune cell types between the PMOP group and control group. (C) Differential expression analysis of immune cell infiltration abundance between PMOP Subtypes 1 and 2. Significance: ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001.

Fig. 8

Correlation analysis of key IR-DEGs with 22 immune cell types. (A) HMOX1 exhibits a significant negative correlation with resting Mast cells. (B) IL1R2 shows a significant positive correlation with activated dendritic cells. (C) SSTR5 demonstrates a significant negative correlation with resting T cells CD4 memory.

Fig. 9

Single-Sample Gene Set Enrichment Analysis (ssGSEA) of the 6 Key IR-DEGs. ssGSEA results for JUN (A), HMOX1 (B), CYSLTR2 (C), TNFSF8 (D), IL1R2 (E) and SSTR5 (F) between high-expression and low-expression groups.

Fig. 10

Diagnostic value of 6 characteristic IR-DEGs in PMOP. (A) ROC curves for each of the 6 characteristic IR-DEGs (JUN, HMOX1, CYSLTR2, TNFSF8, IL1R2, and SSTR5) in the training dataset (GSE56815). (B) ROC curve for the Logistic regression model constructed based on the 6 characteristic IR-DEGs in the GSE56815 dataset. (C) ROC curves for each of the 3 characteristic IR-DEGs in the validation dataset (GSE2208). (D) ROC curve for the Logistic regression model constructed based on the 3 characteristic IR-DEGs in the GSE2208 dataset.

Fig. 11

Potential molecular mechanisms of characteristic IR-DEGs. (A) Construction of the Protein-Protein Interaction (PPI) network based on JUN, HMOX1, TNFSF8, and IL1R2. (B) Generation of the gene-gene interaction network for the key IR-DEGs using the GENEMANIA database. (C) Establishment of the gene-miRNA interaction network based on JUN, HMOX1, TNFSF8, CYSLTR2, and IL1R2. (D) The gene-disease interaction network based on JUN and SSTR5.

Fig. 12

RT-qPCR validation of characteristic IR-DEGs. (A) Analysis of expression differences for three characteristic IR-DEGs (JUN, IL1R2, and CYSLTR2) in the validation set GSE2208. (B) Comparison of microCT results between the OP model group and the control group. (C) Differences in mRNA levels of six characteristic IR-DEGs in peripheral blood between the OP model group and the control group. (D) Differences in mRNA levels of six characteristic IR-DEGs in femoral tissue between the OP model group and the control group. Significance: ns P ≥ 0.05, ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001.