Multi-omics Analysis to Identify Key Immune Genes for Osteoporosis based on Machine Learning and Single-cell Analysis

Abstract

Objective: Osteoporosis is a severe bone disease with a complex pathogenesis involving various immune processes. With the in-depth understanding of bone immune mechanisms, discovering new therapeutic targets is crucial for the prevention and treatment of osteoporosis. This study aims to explore novel bone immune markers related to osteoporosis based on single-cell and transcriptome data, utilizing bioinformatics and machine learning methods, in order to provide novel strategies for the diagnosis and treatment of the disease.

Methods: Single cell and transcriptome data sets were acquired from Gene Expression Omnibus (GEO). The data was then subjected to cell communication analysis, pseudotime analysis, and high dimensional WGCNA (hdWGCNA) analysis to identify key immune cell subpopulations and module genes. Subsequently, ConsensusClusterPlus analysis was performed on the key module genes to identify different diseased subgroups in the osteoporosis (OP) training set samples. The immune characteristics between subgroups were evaluated using Cibersort, EPIC, and MCP counter algorithms. OP's hub genes were screened using 10 machine learning algorithms and 113 algorithm combinations. The relationship between hub genes and immunity and pathways was established by evaluating the immune and pathway scores of the training set samples through the ESTIMATE, MCP-counter, and ssGSEA algorithms. Real-time fluorescence quantitative PCR (RT-qPCR) testing was conducted on serum samples collected from osteoporosis patients and healthy adults.

Results: In OP samples, the proportions of bone marrow-derived mesenchymal stem cells (BM-MSCs) and neutrophils increased significantly by 6.73% (from 24.01% to 30.74%) and 6.36% (from 26.82% to 33.18%), respectively. We found 16 intersection genes and four hub genes (DND1, HIRA, SH3GLB2, and F7). RT-qPCR results showed reduced expression levels of DND1, HIRA, and SH3GLB2 in clinical blood samples of OP patients. Moreover, the four hub genes showed positive correlations with neutrophils (0.65-0.90), immature B cells (0.76-0.92), and endothelial cells (0.79-0.87), while showing negative correlations with myeloid-derived suppressor cells (negative 0.54-0.73), T follicular helper cells (negative 0.71-0.86), and natural killer T cells (negative 0.75-0.85).

Conclusion: Neutrophils play a crucial role in the occurrence and development of osteoporosis. The four hub genes potentially inhibit metabolic activities and trigger inflammation by interacting with other immune cells, thereby significantly contributing to the onset and diagnosis of OP.

Keywords: Bioinformatics; Immunology; Machine learning; Osteoporosis; Single cells analysis.

FIGURE 1

Technology roadmap.

FIGURE 2

Processing of single‐cell data. (A) UMAP plot of GSE147287 cell annotations. (B) Visualization of Top3 genes in various cell types. (C) Violin plots of marker genes for each cell type. (D, E) There are differences in cellular proportions among samples from the high and low bone density groups.

FIGURE 3

hdWGCNA analysis. (A) Annotation of neutrophil subpopulations. (B, C) Distribution of neutrophil subpopulation cell ratios in low and high bone density. (D) The soft threshold selection for hdWGCNA analysis. (E) The yellow module is prominently enriched in clusters 0 and 2.

FIGURE 4

Functional enrichment analysis was performed for the yellow module genes. (A) Gene ontology (GO) enrichment analysis. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis.

FIGURE 5

Cell communication and pseudotime analysis. (A) Cellular communication network diagram. (B) The graph illustrates the receptor interactions of IS_Neutrophils, Other_Neutrophils, and other cells. In the graph, the color represents the magnitude of the binding strength, and larger dots indicate smaller p‐values. (C, D) Pseudotime analysis reveals the temporal milestones of differentiation in different clusters. (E) The yellow module genes are involved in regulating neutrophil differentiation at different time intervals. (F) The transcription factor differences between IS_Neutrophils and Other_Neutrophils are being examined.

FIGURE 6

ConsensusClusterPlus analysis. (A) The training set was divided into cluster 1 (n = 39) and cluster 2 (n = 32) using consensus clustering. (B) The differing expression levels of the yellow module genes in clusters. (C) Volcano plots depicting differential gene expression between clusters. (D) Gene ontology (GO) enrichment analysis of differentially expressed genes. (E) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed genes.

FIGURE 7

Analysis of inter‐subgroup immunomodulatory microenvironments. (A) To analyze the discrepancies in immune characteristics between different clusters, the cibersort, EPIC, and MCPcounter algorithms were utilized. (B) GSVA analysis was performed to evaluate the differences between clusters. (* represents p < 0.05, ** represents p < 0.01, *** represents p < 0.001, **** represents p < 0.0001).

FIGURE 8

Machine learning is used to identify hub genes. (A) Evaluating the performance of validation models based on different combinations of machine learning algorithms. (B) The number of hub genes obtained through different combinations of machine learning algorithms.

FIGURE 9

The relationship between hub genes and immunity and pathways. (A) The correlation between key genes and immune scores was analyzed in the training dataset. (B) The correlation of key genes with pathways in the training set is shown in the figure. The symbols in the figure represent the statistical significance levels: * for p < 0.05, ** for p < 0.01, *** for p < 0.001, **** for p < 0.0001.

FIGURE 10

Analysis of hub gene interactions network. (A) The TF‐hub gene network; (B) The miRNA‐hub gene network.

FIGURE 11

The interaction network analysis of HIRA‐Drug. (A–C) The molecular docking patterns of HIRA with different drugs. The binding free energy of HIRA with Kutizon is −5.41Kcal/mol (A), with clanfenur is −6.64Kcal/mol (B), and with 4‐ethyl‐n‐(3,4,5‐trimethoxybenzyl) aniline, hydrochloride is −3.35Kcal/mol (C). In the figure, purple represents HIRA, cyan represents the drug small molecule, green represents amino acid residues, and red represents the hydrogen bonds between the two.