Abstract Background Postmenopausal osteoporosis is a systemic metabolic disorder typified by an imbalance in bone turnover, where bone resorption supersedes bone formation. This imbalance primarily arises from a decline in bone mass induced by estrogen deficiency, and an elevated risk of fractures resulting from degradation of bone microstructure. Despite recognizing these changes, the precise causative factors and potential molecular pathways remain elusive. In this study, we aimed to identify differentially expressed genes (DEGs), associated pathways, and the role of immune infiltration in osteoporosis, leveraging an integrated bioinformatics approach to shed light on potential underlying molecular mechanisms. Methods We retrieved the expression profiles of [33]GSE230665 from the Gene Expression Omnibus database, comprising 15 femur samples, including 12 postmenopausal osteoporosis samples and 3 normal controls. From the aggregated microarray datasets, we derived differentially expressed genes (DEGs) for further bioinformatics analysis. We used WGCNA, analyzed DEGs, PPI, and conducted GO analysis to identify pivotal genes. We then used the CIBERSORT method to explore the degree of immune cell infiltration within femur specimens affected by postmenopausal osteoporosis. To probe into the relationship between pivotal genes and infiltrating immune cells, we conducted correlation analysis. Results We identified a total of 12,204 DEGs. Among these, 12,157 were up-regulated, and 47 were down-regulated. GO and KEGG pathway analyses indicated that these DEGs predominantly targeted cellular protein localization activity and associated signaling pathways. The protein-protein interaction network highlighted four central hub-genes: RPL31, RPL34, EEF1G, and BPTF. Principal component analysis indicated a positive correlation between the expression of these genes and resting NK cells (as per CIBERSORT). In contrast, the expression of RPL31, RPL34, and EEF1G showed a negative correlation with T cells (gamma delta per CIBERSORT). Conclusions Immune infiltration plays a role in the development of osteoporosis. Keywords: Postmenopausal osteoporosis, WGCNA, Hub genes, Bioinformatics 1. Introduction Osteoporosis is a prevalent condition, especially among postmenopausal women, significantly escalating the risk of fractures. These fractures can lead to severe health complications and even mortality. Chiefly driven by estrogen deficiency, postmenopausal bone loss is the primary factor contributing to the onset of osteoporosis. An estimated 10 % of the global population, including over 30 % of postmenopausal women aged 50 and above, suffer from osteoporosis [[34]1,[35]2]. Thus, understanding the pathological mechanisms of postmenopausal osteoporosis is crucial to developing effective therapeutic strategies. Despite extensive research, the precise mechanism by which E2 loss results in increased bone resorption remains elusive [[36]3]. Other contributing factors include reduced calcium absorption [[37]4,[38]5], a decline in kidney function associated with aging and menopause [[39]6], and impaired vitamin D metabolism [[40]7,[41]8]. Over the past three decades, the understanding of the adaptive immune system's role in postmenopausal osteoporosis development, both in humans and mouse models, has grown significantly. This knowledge of how T-cell-derived cytokines affect bone health led to the birth of osteoimmunology, a term introduced by Arron and Choi in 2000 [[42]9]. Significant strides have been made in understanding the pro-resorptive effects of pro-inflammatory cytokines, especially TNFα and IL-17A produced by T cells. Delving deeper into the complex interplay of immune regulation on bone homeostasis could reveal new, promising targets for therapeutic intervention [[43]10]. Weighted gene co-expression network analysis (WGCNA) has emerged as a powerful tool for processing transcriptomic data. It clusters genes with similar expression patterns into distinct modules and correlates these modules with specific traits, offering a more comprehensive approach than merely examining differentially expressed genes [[44]11]. The adoption of scale-free co-expression network analysis, like WGCNA, has seen a surge across various research domains related to disease [[45][12], [46][13], [47][14]]. With the ongoing advancements in gene detection technologies, such as microarrays and high-throughput sequencing, their potential for identifying key biomarkers crucial for diagnosis, treatment, and prognosis has grown exponentially. Public repositories like the Gene Expression Omnibus (GEO) now serve as indispensable platforms for pinpointing diagnostic and prognostic biomarkers across a myriad of diseases. This study aims to provide foundational insights into the etiology and associated molecular underpinnings of postmenopausal osteoporosis by dissecting differentially expressed genes. Our research seeks to elucidate the role of chromatin regulatory factor-related genes in postmenopausal osteoporosis pathogenesis and to pioneer strategies for potential drug predictions. 2. Materials and methods 2.1. Data sources [48]Fig. 1 outlines the research methodology employed in our study. We sourced the mRNA expression profile microarray [49]GSE230665, a submission by Park et al. ([50]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE230665), from the open-access Gene Expression Omnibus (GEO) database. This dataset encompasses bone samples derived from both postmenopausal osteoporosis-afflicted individuals and healthy controls. The dataset, analyzed using the [51]GPL10332 platform (Agilent-026652 Whole Human Genome Microarray 4 × 44K v2), includes a total of 15 samples: 12 from osteoporosis patients and 3 from control subjects. Fig. 1. [52]Fig. 1 [53]Open in a new tab Workflow chart of the present study. 3. Identification of differentially expressed genes We began with data quality checks and normalization, which included a log transformation to mitigate potential batch effects. The normalization method applied was log2 (X+1). We employed the 'limma' package [[54]15] in R software (Version 4.0.3) to identify differentially expressed genes (DEGs) between postmenopausal osteoporosis patients and controls, setting the criteria at |log FC| ≥2 and a p-value of <0.05. The DEGs were visualized using a volcano plot in R, while the top 30 DEGs were represented in a heatmap generated using the Heatmap package in R. 3.1. Weighted gene Co-expression network analysis We employed WGCNA (as described by Ref. [[55]16] to establish a co-expression network among the DEGs, ensuring adherence to scale-free topology criteria. Starting with an examination of all DEGs using the WGCNA package within R software, we determined the optimal soft thresholding power. With this, a weighted co-expression network was constructed, clustering DEGs into distinct modules, each symbolized by a unique color label. We then explored the correlation between each module and either the postmenopausal osteoporosis patients or the controls. The module with the strongest correlation to postmenopausal osteoporosis patients was identified as the key module and was chosen for further enrichment analysis. 3.2. Gene ontology and pathway enrichment analysis Functional enrichment analysis of the gene set was conducted using the GO annotations [[56]17] from the org.Hs.eg.db R package (version 3.1.0). Mapping and enrichment analyses were performed using the clusterProfiler R package (version 3.14.3). Gene set parameters were set with a minimum of 5 and a maximum of 5000, considering P < 0.05 as statistically significant. Further, the latest KEGG [[57]18,[58]19] gene annotations were procured through its Rest API and underwent a similar enrichment analysis process. Visualization of the derived function and pathway terms was facilitated by Sangerbox [[59]20]. 3.3. Protein–protein interaction network and hub gene identification The PPI network of the crucial module was constructed using the STRING database v11.5 [[60]21] and visualized via Cytoscape software. Significant modules within the PPI network were explored using the MCODE plug-in in Cytoscape, adhering to specified criteria. Hub genes were subsequently identified with the Cytohubba plug-in in Cytoscape. 3.4. Putative signaling pathways involving hub genes and GO analysis 'GeneMania' [[61]22], an expansive web tool, catalogs interaction networks across multiple organisms, assisting in predicting gene functions. Our identified hub genes were fed into 'GeneMania' to construct a presumptive PPI network. All genes from this network were then subjected to GO analysis, with visualization achieved through Sangerbox. 3.5. Assessment of immune cell infiltration The deconvolution algorithm, ‘CIBERSORT’ [[62]23], was employed to analyze gene expression data, discerning the relative proportions of various immune cell types. We conducted an immune infiltration analysis using the CIBERSORT R script and its LM22 gene signature file. This revealed the distribution of 22 immune cell types in postmenopausal osteoporosis samples. Sangerbox visualized these distributions, and Spearman correlation was computed to discern relationships between the immune cell infiltration and the control group, emphasizing specific cell types. 4. Results 4.1. Identifications of DEGs From the comparison between postmenopausal osteoporosis patients and controls, we identified 12,204 DEGs based on criteria: adjusted p < 0.05 and |logFC|≥2.0. Among these, 12,157 genes were upregulated and 47 were downregulated in postmenopausal osteoporosis patients. The volcano plot showcasing these DEGs is provided in [63]Fig. 2A, while a heatmap highlighting the top 30 DEGs is displayed in [64]Fig. 2B. Fig. 2. [65]Fig. 2 [66]Open in a new tab Expression profile of DEGs. A Volcano map of DEGs expression levels. B Heatmap of top 30 DEGs. 4.2. GWCNA analysis Further, the soft-thresholding power β was set to 1, ensuring an optimal average connectivity ([67]Fig. 3A). The DEGs clustered into two modules, namely turquoise and grey, with a set minimum module size of 30. A cluster dendrogram of DEGs is shown in [68]Fig. 3B. Subsequently, we evaluated the correlation of each module with the condition of postmenopausal osteoporosis ([69]Fig. 3C). The turquoise module emerged as the most significantly correlated with the condition (correlation: 0.72, p < 0.0001), housing 133 DEGs (MM = 0.9, GS = 0.8). Fig. 3. [70]Fig. 3 [71]Open in a new tab WGCNA of DEGs. A. Estimation of the soft thresholding value for a scale-free co-expression network. B.Cluster dendrogram of all DEGs. C. Correlation between each module and postmenopausal osteoporosis patients or controls. 4.3. Functional enrichment analysis The turquoise module's 133 DEGs underwent functional enrichment analysis in R. The GO analysis revealed significant enrichment in 16 biological processes (BP), 73 cellular components (CC), and 77 molecular functions (MF). The top enrichments for BP, CC, and MF are delineated in [72]Fig. 4A–C. Prominent categories like 'cellular protein localization,' 'nuclear lumen,' and 'structural molecule activity' were discerned. Further, KEGG analysis linked the turquoise module to two pathways: 'Ribosome' and 'Chronic myeloid leukemia' ([73]Fig. 4D). Fig. 4. [74]Fig. 4 [75]Open in a new tab Gene ontology and KEGG enrichment analysis. A. Biological process. B. Cellular component. C. Molecular function. D. KEGG enrichment analysis. 4.4. PPI network Construction, Modular analysis, and hub gene analysis The gene interactions within the turquoise module were examined, forming a protein–protein interaction network via the STRING database. This network consisted of 120 nodes and 75 edges. The network's most significant module is depicted in [76]Fig. 5A. Using the MCODE plugin in Cytoscape, three modules were identified. Module 1 (score: 6.667) had 7 nodes and 20 edges ([77]Fig. 5B); Module 2 (score: 5.000) featured 5 nodes and 10 edges ([78]Fig. 5C); Module 3 (score: 3.000) contained 3 nodes and 3 edges ([79]Fig. 5D). [80]Table 1 enumerates the top ten pivotal genes, as deduced by five algorithms (MCC, DMNC, MNC, Degree, and EPC) in the cellular mapping plugin. Overlapping hub genes across these algorithms, namely RPL31, RPL34, EEF1G, and BPTF, were corroborated using a Venn diagram ([81]Fig. 5E). Fig. 5. [82]Fig. 5 [83]Open in a new tab PPI network and hub gene. A. PPI network. B, C, D. The three most significant modules. E. The overlapped hub genes from different algorithms. Table 1. Top ten hub genes obtained by five algorithms of Cytohubba. MNC MCC EPC DMNC GEGREE RPL34 RPL5 CREBBP RPL34 CREBBP RPL36 RPL31 BPTF RPL36 BPTF RPL5 RPL34 EED RPL5 RPL5 RPS19 RPL36 KMT2C RPS19 RPL31 RPL31 RPS19 SETD1A RPL31 EEF1G EEF1G EEF1G MED13L EEF1G RPL34 BPTF SRP14 MECOM SRP14 RPL36 SRP14 BPTF RPL31 MED13L RPS19 KMT2C KMT2C RPL34 CREBBP KMT2C [84]Open in a new tab Construction and Analysis of Putative RPL31, RPL34, EEF1G, and BPTF Protein–Protein Interaction (PPI) Network. By leveraging the 'GeneMania' tool, we built a putative protein–protein interaction network comprising 24 genes, which notably included the hub genes: RPL31, RPL34, EEF1G, and BPTF. Illustrated in [85]Fig. 6A, this PPI network encompasses a total of 2364 links. Subsequent gene ontology analysis revealed that these 24 genes predominantly clustered within 99 biological processes, 50 cellular components, and 13 molecular functions. These enrichments are visualized in [86]Fig. 6B–D. It is inferred that the hub genes RPL31, RPL34, EEF1G, and BPTF, which are associated with postmenopausal osteoporosis, play roles in 'translation', 'cytosol', and 'nucleic acid binding'. This suggests the possibility of significant nucleic acid alterations in the pathogenesis of postmenopausal osteoporosis.3.6 Evaluation of immune cell infiltration and immune cell correlation analysis. Fig. 6. [87]Fig. 6 [88]Open in a new tab Putative RPL31,RPL34, and BPTF protein–protein interaction network and GO analysis A. PPI network. B. Biological process. C. Cellular component. D. Molecular function. 4.5. Estimation of leukocyte Subpopulation proportions and their association with hub genes To decipher the immune landscape associated with postmenopausal osteoporosis, we employed the CIBERSORT tool. This allowed us to estimate the proportions of 22 distinct leukocyte subtypes within femur samples. These subtypes spanned a range from native B cells and different T cell subtypes to various macrophages and granulocytes ([89]Fig. 7A).Upon visualization using violin plots, a striking differential immune cell infiltration pattern emerged between postmenopausal osteoporosis patients and control subjects. Specifically, heightened levels of T_cells_gamma_delta_CIBERSORT and NK_cells_resting_CIBERSORT were evident in osteoporosis-afflicted samples, indicating their potential role in disease pathophysiology ([90]Fig. 7B).Furthermore, we sought to understand the interplay between these altered immune cell populations and our identified hub genes. Spearman's correlation analysis revealed an intriguing association. NK_cells_resting_CIBERSORT exhibited a positive correlation with the expression levels of hub genes RPL31, RPL34, EEF1G, and BPTF. Contrastingly, T_cells_gamma_delta_CIBERSORT was inversely related to the expression of RPL31, RPL34, and EEF1G, indicating their potential antagonistic roles in the context of postmenopausal osteoporosis ([91]Fig. 7C). Fig. 7. [92]Fig. 7 [93]Open in a new tab Immune cell infiltration analysis. A. Heat map of relative proportions of 22 infiltrated immune cells in patients with Postmenopausal osteoporosis. B. Violin chart of the abundance of each type of immune cell infiltration in Postmenopausal osteoporosis and control groups. C. The correlation analysis of hub-genes and NK_cells_resting_CIBERSORT, T_cells_gamma_delta_CIBERSORT in Postmenopausal osteoporosis. 5. Discussion Osteoporosis, particularly evident after menopause, is significantly driven by a decline in estrogen levels. The risk of developing postmenopausal osteoporosis is exacerbated by several factors: aging, genetic factors, smoking, low body weight, and certain diseases or medications that compromise bone health. When treatment becomes necessary, various approved remedies are available, including estrogen agonists/antagonists, bisphosphonates, RANK ligand inhibitors, parathyroid hormone receptor agonists, and sclerostin inhibitors. The World Health Organization's report highlights the severity of the issue, noting that postmenopausal women face at least a 30 % risk, if not closer to 40 %, of suffering osteoporotic fractures. This stark statistic underscores the profound health implications of osteoporotic bone injuries for middle-aged and elderly women, emphasizing the urgency of addressing this significant clinical challenge. In our study, we employed the WGCNA method on the mRNA expression profile [94]GSE230665 obtained from the GEO database. Utilizing this cutting-edge methodology, we categorized all 12,204 differentially expressed genes, as recognized by the 'limma' algorithm, into two unique modules. Notably, among these, 133 genes within the turquoise module displayed a pronounced correlation with postmenopausal osteoporosis, evidenced by a correlation score of 0.72 and a p-value less than 0.0001. Further analysis revealed that these genes were predominantly associated with 'cellular protein localization,' 'nuclear lumen,' and 'structural molecule activity.' A deeper dive into the Protein-Protein Interaction (PPI) network identified four pivotal genes: RPL31, RPL34, EEF1G, and BPTF. Functional enrichment analysis of these genes spotlighted their significant involvement in processes like 'translation,' 'cytosol,' and 'nucleic acid binding'. Ribosomal Proteins (RPs) serve as the foundational components of the ribosome complex. They play pivotal roles in cell growth and proliferation by assembling systematically. Not just vital for optimizing protein synthesis, these RPs also function as surveillance mechanisms for cellular health [[95]24]. A dataset on rheumatoid arthritis unveiled 113 differentially expressed genes (DEGs), with 104 being up-regulated and 9 down-regulated. Notably, seven pivotal genes were identified from this set, including RPL34, as revealed through PPI network analysis [[96]25]. Delving deeper into specific RPs: RPL31 (also termed L31): Situated on human chromosome 2 (2q11.2), it's integral to the 60S large ribosomal subunit. RPL31 partakes in numerous processes, ranging from ribosome self-assembly and protein synthesis to DNA repair and tumorigenesis [[97]26,[98]27,[99]28]. A study revealed that MiR-483-5p targets RPL31, consequently inhibiting the RAS/MEK/ERK signaling pathway, and thereby affecting osteogenic differentiation. This positions RPL31 as a promising therapeutic target for osteoporosis and bone defects [[100]29]. RPL34: A member of the L34E RP family, it is located on the long arm of human chromosome 4. Its key function in forming the 60S subunit links it to the progression of various malignant tumors [[101]30]. EEF1: This Eukaryotic Elongation Factor complex is composed of six subunits. Its aberrant regulation can lead to alterations in the epigenomic landscape, influencing the overall gene transcription profile and contributing to cancer onset. These unique molecular features offer potential pathways for targeted drug therapies and personalized medicine approaches [[102]31] iterge-Sut,2019. BPTF: As the largest component of the chromatin remodeling complex, BPTF plays a pivotal role in regulating DNA accessibility and gene expression [[103][32], [104][33], [105][34]]. Its significance is further underscored by its identification as a potential biomarker gene for rheumatoid arthritis [[106]35]. Estrogen deficiency-induced bone loss is a multifaceted phenomenon, with the immune system playing a more significant role than the direct effects of estrogen on bone cells [[107]36]. Recent clinical findings highlight this intricate interplay. A cross-sectional study showcased higher white blood cell counts in postmenopausal women than their premenopausal peers, pointing towards an increase in overall lymphocyte and monocyte counts [[108]37]. This observation is in line with prior experiments showing that estrogen deficiency boosts T cell outputs from the thymus into the peripheral blood [[109]38]. Furthermore, in postmenopausal osteoporosis, there's a peculiar pattern: while B-lymphocyte levels in the bone marrow are normal or even reduced, there's a spike in individual G-CSF secretion. Dive deeper into this phenomenon, studies involving ovariectomized (OVX) mice have shown that, without estrogen, B cells tend to proliferate more. This proliferation seems to be spurred by an increase in CXCL12 signaling. These cells also ramp up their G-CSF secretion [[110]39,[111]40]. Interestingly, a specific strain of mast cell-deficient mice (Mcpt5-Cre R-DTA flox) appears to be resistant to bone loss usually seen after ovariectomy. This suggests that, in the absence of estrogen, mast cells may promote osteoclastogenesis, leading to bone loss [[112]41]. Our own study adds another layer to this understanding. We've found an augmented presence of T cells and NK cells during the progression of postmenopausal osteoporosis. While T cells showed a negative association with three of our identified core genes, NK cells positively correlated with all four. Notably, the exact role of NK cells in postmenopausal osteoporosis remains a relatively uncharted territory, underscoring the need for more comprehensive research in this domain. While this study provides important insights, we recognize its limitations. Our primary reliance on bioinformatics analysis of public datasets may not fully capture real-world physiological contexts. The use of just one dataset for postmenopausal osteoporosis potentially narrows the breadth of our findings. To address these concerns, we intend to monitor database updates for potential enrichments to our data and to obtain relevant clinical samples for further gene sequencing. Importantly, the lack of experimental validation leaves our analysis not entirely conclusive. Moving forward, we aim to validate and expand upon our findings using animal models and in-depth cellular experimentation. 6. Conclusion This study has identified four pivotal hub genes believed to be instrumental in the pathophysiology of postmenopausal osteoporosis. Notably, we observed an elevation in specific immune cells, particularly T cells and NK cells, in affected patients. These immune cells displayed notable interactions with the aforementioned hub genes. The interplay between these elements presents an intricate and potentially significant mechanism, warranting further in-depth research to conclusively elucidate their roles and interactions in the context of postmenopausal osteoporosis. Funding statement This work was supported by Tianjin Key Medical Discipline(Specialty) Construction Project (TJYXZDXK-026A) and Hebei Province High Level Talent Funding Project (C20221122). Data availability statement Data will be made available on request. No additional information is available for this paper. CRediT authorship contribution statement Xiaoli Zhou: Writing – original draft. Yang Chen: Writing – review & editing. Zepei Zhang: Data curation. Jun Miao: Visualization. Guangdong Chen: Formal analysis. Zhiyong Qian: Validation. Declaration of competing interest The authors declare no conflict of interest. Contributor Information Jun Miao, Email: mj6688@tju.edu.cn. Guangdong Chen, Email: victordongdong@126.com. Abbreviations list BP Biological processes CC Cellular components DEGs Differentially expressed genes EEF1 Eukaryotic Elongation Factor complex 1 GEO Gene Expression Omnibus GO Gene Ontology KEGG Kyoto Encyclopedia of Genes and Genomes MCODE Molecular Complex Detection MF Molecular functions PPI Protein-Protein Interaction RP Ribosomal Proteins STRING Retrieval Interacting Genes WGCNA Weighted Gene Co-Expression Network Analysis References