Abstract

   Steroid-induced osteonecrosis of the femoral head (SONFH) is a
   debilitating condition caused by long-term corticosteroid use, leading
   to impaired blood flow and bone cell death. The disruption of cellular
   processes and promotion of apoptosis by endoplasmic reticulum stress
   (ERS) is implicated in the pathogenesis of SONFH. We identified
   ERS-associated genes in SONFH and investigated their potential as
   therapeutic targets. We analysed the [40]GSE123568 GEO dataset to
   identify differentially expressed genes (DEGs) related to ERS in SONFH.
   We conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and
   Genomes (KEGG) pathway enrichment analyses, identified hub genes by
   protein-protein interaction (PPI) analyses, and evaluated their
   functions by gene set enrichment analysis (GSEA). We constructed
   mRNA–miRNA networks, identified potential therapeutics, and assessed
   immune cell infiltration. We performed cross-validation using the GEO
   dataset [41]GSE74089, qRT-PCR on clinical samples from patients with
   SONFH and controls, and a receiver operating characteristic (ROC) curve
   analysis to assess the diagnostic performance of the hub genes. We
   identified 195 ERS-related genes in SONFH, which were primarily
   involved in oxidative stress, immune responses, and metabolic pathways.
   The PPI network suggested CXCL8, STAT3, IL1B, TLR4, PTGS2, TLR2, CASP1,
   CYBB, CAT, and HOMX1 to be key hub genes, which were shown by GSEA to
   be involved in biological pathways related to metabolism, immune
   modulation, and cellular integrity. We also identified 261 microRNAs
   (miRNAs) as well as drugs such as dibenziodolium and
   N-acetyl-L-cysteine that modulated inflammatory responses in SONFH.
   Twenty-two immune cell subtypes showed significant correlations, such
   as a positive correlation between activated mast cells and Tregs, and
   patients with SONFH had fewer dendritic cells than controls. The hub
   genes CYBB and TLR4 showed significant correlations with M1 macrophages
   and CD8 T cells, respectively. Cross-validation and qRT-PCR confirmed
   the upregulation of STAT3, IL1B, TLR2, and CASP1 in patients with
   SONFH, validating the bioinformatics findings. An ROC curve analysis
   confirmed the diagnostic potential of the hub genes. The top 10 hub
   genes show promise as ERS-related diagnostic biomarkers for SONFH. We
   discovered that 261 miRNAs, including hsa-miR-23, influence these genes
   and identified potential therapeutics such as dibenziodolium and
   simvastatin. Immune profiling indicated altered immune functions in
   SONFH, with significant correlations among immune cell types.
   Validation confirmed the upregulation of STAT3, IL1B, TLR2 and CASP1,
   which had diagnostic potential. The findings suggest potential
   diagnostic markers and therapeutic targets for SONFH.

   Subject terms: Computational biology and bioinformatics, Biomarkers,
   Musculoskeletal system, Osteoimmunology

Introduction

   Osteonecrosis of the femoral head (ONFH) is a debilitating condition
   characterised by ischaemia and cellular death in the femoral head,
   leading to structural and functional changes in the bone^[42]1–[43]3.
   The aetiology of OFNH can be traumatic or non-traumatic^[44]4. Hormone
   use, chronic alcohol consumption, and a history of trauma are
   implicated as triggers of ONFH^[45]5. Steroid use is the most prevalent
   non-traumatic cause of ONFH^[46]6.

   Steroid-induced ONFH (SONFH) occurs predominantly as a result of
   excessive use of corticosteroids in patients with immune-related
   diseases^[47]7,[48]8. The effects of glucocorticoids (GCs) are
   important in SONFH. GCs can induce lipid metabolism disorders, leading
   to intraosseous hypertension, impairing the blood supply to the femoral
   head and initiating necrosis^[49]9. In addition, GCs can directly
   induce apoptosis of osteoblasts and osteocytes, which are responsible
   for bone formation and maintenance, thus accelerating bone loss and
   necrosis^[50]10. Moreover, GCs influence autophagy, which is essential
   for maintaining cellular homeostasis, by altering the balance between
   autophagy activation and inhibition^[51]11. This may contribute to the
   impairment of osteocyte function and increase susceptibility to ONFH.
   Furthermore, GCs reportedly regulate ferroptosis, a type of cell death
   driven by iron-dependent lipid peroxidation. Excess GC exposure
   promotes ferroptosis, which may aggravate bone loss and contribute to
   the progression of SONFH^[52]12. Critically, GCs activate endoplasmic
   reticulum stress (ERS), exacerbating the disruption of osteogenesis and
   angiogenesis, which are key factors in disease development^[53]13.

   ERS occurs when the ER (the site of protein folding, lipid synthesis,
   and calcium storage) is overwhelmed by unfolded or misfolded proteins,
   triggering the unfolded protein response (UPR)^[54]13. The UPR is
   mediated by three key pathways: IRE1, which enhances the ability of the
   ER to manage misfolded proteins by splicing XBP1 mRNA; ATF6, which upon
   ER stress travels to the Golgi apparatus to be cleaved and induce the
   expression of UPR target genes; and PERK, which phosphorylates eIF2α to
   reduce protein synthesis and lessen the ER load, and activates stress
   response genes^[55]14,[56]15. These pathways restore ER function by
   reducing protein misfolding, enhancing protein folding, and degrading
   misfolded proteins. If this balance is not restored, prolonged ERS can
   lead to apoptosis, exacerbating bone cell dysfunction and accelerating
   the progression of SONFH^[57]16.

   In this study, we used bioinformatics methods to identify
   differentially expressed genes (DEGs) associated with ERS in SONFH. We
   analysed the roles of key hub genes, investigated their interactions
   with regulatory microRNAs (miRNAs), and identified therapeutics that
   could target these genes. In addition, we examined immune cell
   infiltration to understand the immunological landscape of SONFH. We
   validated the findings using external datasets and clinical samples.
   Our results will facilitate the development of targeted therapeutics
   for SONFH.

Materials and methods

Bioinformatics

Data collection

   We downloaded the [58]GSE123568 dataset from the Gene Expression
   Omnibus (GEO) database ([59]https://www.ncbi.nlm.nih.gov/geo/), and
   compared gene expression levels in 30 patients with SONFH and 10
   non-affected controls post-steroid treatment^[60]17. Next, we extracted
   genes associated with ERS from the GeneCards database
   ([61]https://www.genecards.org)^[62]18. Using a relevance threshold
   score of 10, we identified 2176 ERS-related genes (ERGs).

Analysis of differentially expressed genes

   We used DESeq2 in R to identify genes in the [63]GSE123568 dataset with
   significant changes in expression, defined as an absolute
   fold-change > 0.75 and an adjusted p < 0.05^[64]19; these were termed
   DEGs. The DEGs intersecting with 2176 ERGs were classified as
   differentially expressed ERS-related genes (DE-ERGs).

Functional enrichment analysis

   We conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and
   Genomes (KEGG) analyses to evaluate the functions of the
   DE-ERGs^[65]20–[66]22. GO analysis categorises DE-ERGs into Biological
   Processes (BP), Cellular Components (CC), and Molecular Functions (MF).
   KEGG analysis maps DE-ERGs to broader biological pathways. We used the
   clusterProfiler package to identify significant enrichments (adjusted
   p < 0.05 and Q < 0.2) and visualised the results with ggplot2^[67]23.

Screening of hub genes

   We used the STRING database ([68]https://string-db.org/) to analyse
   protein–protein interactions (PPIs) related to DE-ERGs^[69]24. The
   dataset was imported into Cytoscape, and we used the cytoHubba plugin
   to identify hub genes with crucial connections in the ERS PPI
   network^[70]25. We used the MCODE plugin in Cytoscape to identify
   significant clusters and refined the focus to the top 10 hub genes.

Single-gene gene set enrichment analysis

   We conducted a gene set enrichment analysis (GSEA) to evaluate the
   associations between the hub genes with biological pathways or
   processes. Using the Molecular Signatures Database (MSigDB)
   ([71]https://www.gseamsigdb.org/gsea/msigdb/index.jsp), we linked the
   hub genes to biological functions^[72]26. This approach is underpinned
   by stringent criteria (adjusted p < 0.05 and false discovery rate
   [FDR] < 0.25). We used the ggplot2 package in R to visualise the
   contributions of the hub genes to ERS and other stress-response
   mechanisms.

Construction of the mRNA–miRNA network

   We identified miRNAs that may regulate the top 10 genes implicated in
   ERS using miRWalk ([73]http://mirwalk.umm.uni-heidelberg.de/)^[74]27.
   By visualising these interactions in Cytoscape, we constructed an
   mRNA–miRNA regulatory network.

Targeted drug prediction

   We used Enrichr ([75]https://maayanlab.cloud/Enrichr) to identify drugs
   linked to key genes in ERS^[76]28,[77]29. We input the top 10 hub genes
   into Enrichr, applied the Drug Signatures Database (DSigDB), and
   prioritised the results by p-values to identify the top 10 drugs
   associated with the hub genes.

Analysis of immune cell infiltration

   We used the CIBERSORT algorithm ([78]https://cibersortx.stanford.edu/)
   to analyse immune cell infiltration^[79]30. Using a permutation count
   of 100 and a p < 0.05, we evaluated the populations of immune cell
   types. We used the vioplot package to visualise the abundances of the
   immune cell types. The corrplot package was employed to generate
   correlation matrices to provide insight into the relationships between
   hub genes and immune cell types.

Cross-validation using an external dataset

   To validate the findings from the [80]GSE123568 dataset, we performed
   cross-validation^[81]31 of DE-ERG expression levels in the [82]GSE74089
   dataset. This dataset includes gene expression profiles of hip
   cartilage from patients with necrosis of the femoral head. By comparing
   the expression levels of the DE-ERGs in [83]GSE74089 and [84]GSE123568,
   we validated the identified hub genes. A value of p < 0.05 was
   considered indicative of statistical significance.

Validation of clinical samples

Clinical samples

   This study complied with the Declaration of Helsinki and was approved
   by the Medical Ethics Committee of Qingdao University Affiliated
   Hospital. All of the participants provided informed consent. The
   inclusion criteria were adults 18–60 years of age diagnosed via
   radiographic evidence with femoral neck fractures or SONFH. The
   exclusion criteria were previous conditions that could influence immune
   response or gene expression, such as systemic infections or autoimmune
   diseases; immunosuppressive therapy within 6 months; or a history of
   metabolic bone diseases unrelated to our condition of interest, e.g.,
   osteoporosis. We collected peripheral blood samples from three patients
   with SONFH and three healthy adults (following steroid administration).

Quantitative real-time polymerase chain reaction (qRT-PCR)

   Peripheral blood was collected into EDTA-treated tubes and centrifuged
   to segregate the cellular fraction from plasma. White blood cells were
   lysed using TRIzol reagent, releasing RNA and denaturing proteins and
   DNA. RNA was isolated from the aqueous phase, precipitated using
   alcohol, and the pellet was washed with ethanol to remove impurities.
   After resuspension in RNase-free water, RNA purity and concentration
   were determined. We used the qPCR HiScript III RT SuperMix Kit to
   convert RNA into cDNA according to the manufacturer’s protocol. cDNA
   was used for qRT-PCR with the ChamQ Universal SYBR qPCR Master Mix Kit,
   with ACTB mRNA as the normalisation control. Each sample underwent
   qRT-PCR in triplicate, together with three biological replicates per
   group. The primer sequences are listed in Table [85]1.

Table 1.

   Primers used for real-time PCR.
   Gene  Species Direction Sequence (5ʹ to 3ʹ)
   CXCL8 Human   Forward   ACTTTCAGAGACAGCAGAGCACAC
                 Reverse   CACACAGTGAGATGGTTCCTTCCG
   STAT3 Human   Forward   GAGGCAGGAGAATCGCTTGAACC
                 Reverse   TCTCAGACTGTCGCCCAGGATG
   IL1B  Human   Forward   GACCTGGACCTCTGCCCTCTG
                 Reverse   GCCTGCCTGAAGCCCTTGC
   TLR4  Human   Forward   CATCATTGGTGTGTCGGTCCTCAG
                 Reverse   AGCCAGCAAGAAGCATCAGGTG
   PTGS2 Human   Forward   ACAGCCAGACGCCCTCAGAC
                 Reverse   ATGGGTGGGAACAGCAAGGATTTG
   TLR2  Human   Forward   GTGACCAAGTGAAGGCAGGAAGAC
                 Reverse   AGCAGGCAGGAGGGGTGTTG
   CASP1 Human   Forward   CACACCGCCCAGAGCACAAG
                 Reverse   TCCCACAAATGCCTTCCCGAATAC
   CYBB  Human   Forward   TTCCAGTGCGTGCTGCTCAAC
                 Reverse   ACTCGGGCATTCACACACCATTC
   CAT   Human   Forward   CGGCAAGGGAGAAGGCAAATCTG
                 Reverse   GATGAGCGGGTTACACGGATGAAC
   HMOX1 Human   Forward   TGCCAGTGCCACCAAGTTCAAG
                 Reverse   TGTTGAGCAGGAACGCAGTCTTG
   [86]Open in a new tab

Validation of hub genes by receiver operating characteristic analysis

   Using the R package pROC, we conducted ROC curve analyses on the
   [87]GSE123568 dataset to evaluate the diagnostic accuracy of hub genes
   for SONFH^[88]32. We assessed model performance by calculating areas
   under the curve (AUCs); AUCs > 0.6 indicated statistical significance.

Results

DE-ERGs in SONFH

   We identified 1220 DEGs, including 754 upregulated and 466
   downregulated genes, in the [89]GSE123568 dataset (Fig. [90]1A). By
   intersecting these DEGs with ERS-related genes scoring > 10 in the
   GeneCards database, we identified 195 DE-ERGs (Fig. [91]1B). The
   expression levels of the top 100 intersecting genes of the SONFH and
   control groups are shown in Fig. [92]1C.

Fig. 1.

   [93]Fig. 1
   [94]Open in a new tab

   Analysis and identification of differentially expressed endoplasmic
   reticulum stress-related genes (DE-ERGs) in SONFH and control samples.
   (A) Volcano map of the distribution of the 1220 DEGs in the
   [95]GSE123568 dataset. Red, upregulated DEGs; blue, downregulated DEGs.
   (B) Venn diagram of 195 DE-ERGs. (C) Heatmap of changes in expression
   of the top 100 DE-ERGs. Blue, downregulation; red, upregulation.

GO and KEGG enrichment analyses of DE-ERGs

   The GO and KEGG analyses identified 1987 GO terms and 117 KEGG pathways
   related to SONFH. The GO analysis yielded 1682 BPs, 137 CCs, and 158
   MFs (Fig. [96]2), including BPs such as response to oxidative stress
   (OS), peptide-hormone signalling, and regulation of defence responses.
   CCs included structural features such as membrane microdomains and
   mitochondrial outer membranes. MFs were dominated by kinase activity
   and binding mechanisms. KEGG pathway analysis linked these DE-ERGs to
   the lipid and atherosclerosis, NOD-like receptor signalling, and FoxO
   signalling pathways.

Fig. 2.

   [97]Fig. 2
   [98]Open in a new tab

   Enrichment analyses of DE-ERGs. (A,B) DE-ERGs related to Biological
   Process (BP), Cellular Component (CC), and Molecular Function (MF)
   according to GO analysis. (C,D) KEGG pathway enrichment analysis.

Construction of the PPI network and identification of hub genes

   We constructed a PPI network for the 195 DE-ERGs using the STRING
   database with visualisation in Cytoscape (Fig. [99]3). The MCODE plugin
   revealed three critical gene clusters. Cluster 1 included genes
   involved in inflammation and immune response, Cluster 2 contained genes
   linked to OS and innate immunity, and Cluster 3 included lipid
   metabolism enzymes. The cytoHubba plugin identified CXCL8, STAT3, IL1B,
   TLR4, PTGS2, TLR2, CASP1, CYBB, CAT, and HMOX1 as hub genes. Nine genes
   were upregulated, and one gene (CAT) was downregulated.

Fig. 3.

   [100]Fig. 3
   [101]Open in a new tab

   Protein–protein interaction (PPI) network and top 10 hub genes. (A) PPI
   network of 195 DE-ERGs. (B) The MCODE plug-in identified three
   important clusters. (C) Identification of the top 10 hub genes using
   the CytoHubba plug-in. (D) Up- and downregulated genes among the top 10
   hub genes.

   In the upregulated gene network, coexpression interactions predominated
   at 66.92%, followed by colocalisation at 18.35%, and physical
   interactions at 15.43%. Predicted interactions were noted at 4.46%,
   with pathway involvement at 4.14%, genetic interactions at 0.46%, and
   shared protein domains at 0.23%. The network of the downregulated gene,
   CAT, was characterised by physical interactions, which constituted
   77.64% of the connections. Coexpression accounted for 8.1%, with
   predicted interactions at 5.37%, colocalisation at 3.63%, genetic
   interactions at 2.87%, pathway involvement at 1.88%, and shared protein
   domains at 0.60%.

Potential mechanisms of hub genes explored by GSEA

   We performed GSEA for each of the hub genes (Fig. [102]4). CXCL8 showed
   enrichment in Thyroid Cancer and Galactose Metabolism, and negative
   enrichment in SNARE Interactions in Vesicular Transport and Proteasome.
   STAT3 showed negative enrichment in Epithelial Cell Signalling in
   Helicobacter pylori Infection and Prostate Cancer. IL1B showed positive
   enrichment in Insulin Signalling Pathway and NOD-like Receptor
   Signalling Pathway, and TLR4 exhibited negative enrichment in T Cell
   Receptor Signalling Pathway and ERBB Signalling Pathway. GSEA of PTGS2
   indicated negative enrichment in, for instance, NOD-like Receptor
   Signalling Pathway and Apoptosis. TLR2 was negatively enriched in
   Vibrio cholerae Infection and Toll-like Receptor Signalling Pathway.
   CASP1 showed a mix of positive enrichment in pathways such as Pentose
   and Glucuronate Interconversions and negative enrichment in Lysosome
   and SNARE Interactions in Vesicular Transport. GSEA of CYBB revealed
   negative enrichment in Autoimmune Thyroid Disease and Lysosome. CAT was
   positively enriched in Valine, Leucine, and Isoleucine Degradation and
   Glyoxylate and Dicarboxylate Metabolism. HMOX1 exhibited negative
   enrichment in Selenoamino Acid Metabolism and Glycerolipid Metabolism.

Fig. 4.

   [103]Fig. 4
   [104]Open in a new tab

   Gene set enrichment analysis (GSEA) of the top 10 hub genes. Colours
   above the image are labelled with the names of the enrichment entries
   corresponding to the colours in the image.

Hub gene-miRNA regulatory network and drugs

   The mRNA–miRNA network showed that 261 miRNAs modulated the expression
   of the top 10 hub genes associated with ERS (Fig. [105]5). These
   included hsa-miR-23, which is reportedly correlated with osteogenic and
   adipogenic functions^[106]33–[107]36. In addition, dibenziodolium and
   N-acetyl-L-cysteine, together with simvastatin and curcumin (which
   modulated inflammatory responses), potentially interacted with the hub
   genes.

Fig. 5.

   [108]Fig. 5
   [109]Open in a new tab

   Hub gene-miRNA regulatory network and targeted drugs. (A) Interaction
   network between hub genes and targeted miRNAs. (B) Drugs identified in
   the analysis of the top 10 hub genes.

Immune infiltration and hub gene-immune cell interactions

   Immune infiltration analysis using the [110]GSE123568 dataset suggested
   the presence of diverse immune-cell subpopulations. The relative
   percentages of 22 immune-cell subtypes are shown in Fig. [111]6. The
   correlation analysis indicated significant interactions between
   immune-cell types. For instance, activated mast cells were strongly
   positively correlated with Tregs and resting dendritic cells and
   resting NK cells were negatively correlated with gamma-delta T cells.

Fig. 6.

   [112]Fig. 6
   [113]Open in a new tab

   Immune infiltration and hub gene-immune cell interactions. (A) Relative
   percentages of 22 subpopulations of immune cells in 40 samples from the
   [114]GSE123568 dataset. (B) Spearman correlation analysis among immune
   cells. Red, positive correlation; blue, negative correlation. (C)
   Boxplots of the expression of immune cells between the two groups. (D)
   Correlation map of hub genes and immune infiltration. Yellow, positive
   correlation; blue, negative correlation. The darker the colour the
   stronger the correlation (*p < 0.05).

   Immune-landscape analysis showed that the proportion of dendritic cells
   was significantly smaller in the SONFH group compared to the control
   group, suggesting suppression or deficit in antigen presentation or
   immune priming in SONFH.

   We analysed the correlation between the hub genes identified in the PPI
   network and immune cell infiltration. The gene CYBB showed a
   significant positive correlation with M1 macrophages and a negative
   correlation with activated NK cells, indicating an inverse relationship
   with adaptive immune responses. TLR4 showed a significant negative
   correlation with CD8 T cells.

Validation of hub genes in the [115]GSE123568 dataset

   The cross-validation results confirm that the expression levels of the
   10 hub genes in the [116]GSE74089 dataset align with those in the
   [117]GSE123568 dataset. Relative to the control group, the SONFH group
   exhibited upregulation of CXCL8, STAT3, IL1B, TLR4, PTGS2, TLR2, CASP1,
   CYBB, and HMOX1, and downregulation of CAT (Fig. [118]7). qRT-PCR
   analysis of samples from three patients with SONFH and three healthy
   adults confirmed these findings at the mRNA level, reinforcing the
   reproducibility of the [119]GSE123568 data. Compared to the control
   group, there was significant upregulation of STAT3 (p = 0.0002), IL1B
   (p = 0.0396), TLR2 (p = 0.0054), and CASP1 (p = 0.0264) in the SONFH
   group, which agrees with the bioinformatics analysis (Fig. [120]8). To
   validate the predictive significance of the 10 hub genes, we performed
   an ROC curve analysis on the [121]GSE123568 dataset. All of the hub
   genes had statistically significant AUC values > 0.6 (Fig. [122]9).
   PTGS2 (AUC = 0.805) showed the greatest diagnostic efficacy, followed
   by STAT3 (AUC = 0.815).

Fig. 7.

   [123]Fig. 7
   [124]Open in a new tab

   Violin plot of 10 hub genes in SONFH and control samples in the
   [125]GSE74089 dataset.

Fig. 8.

   [126]Fig. 8
   [127]Open in a new tab

   Validation by qRT-PCR of the top 10 hub genes in SONFH and control
   samples. Experiments were performed in triplicate; results are
   means ± SD (*p < 0.05).

Fig. 9.

   [128]Fig. 9
   [129]Open in a new tab

   Diagnostic receiver operating characteristic (ROC) curves of 10 hub
   genes in SONFH samples and control samples.

Discussion

   Our findings provide insight into the molecular framework of SONFH by
   identifying and characterising DE-ERGs. Analysis of the [130]GSE123568
   dataset revealed a large number of DE-ERGs, thereby enhancing our
   understanding of the cellular mechanisms triggered by steroid exposure.
   The study that originally generated the [131]GSE123568 dataset focused
   on identifying dynamic biomarkers for different stages of non-traumatic
   ONFH (NONFH) using a transcriptional regulatory network^[132]17.
   Specifically, it highlighted the alterations in gene expression across
   various stages of the disease and identified key regulatory genes and
   pathways involved in NONFH progression. However, our study specifically
   investigates the role of ERS in SONFH, which was not addressed in the
   original study. Our research fills this gap by specifically
   investigating the intersection of ERS-related genes with the DEGs
   identified from the [133]GSE123568 dataset. The identification of 1220
   DEGs (754 upregulated and 466 downregulated) highlights the effects of
   steroids on gene expression in SONFH. Among these, the intersection
   with ERS-related genes resulting in 195 DE-ERGs suggests a significant
   subset of genes that may play crucial roles in the pathophysiology of
   SONFH through ERS mechanisms. This focused analysis provides a more
   detailed understanding of how ERS contributes to the cellular and
   molecular disruptions observed in SONFH, thereby offering new insights
   that were not explored in the original study.

   GO and KEGG enrichment analyses provided insight into the BP and
   pathways enriched for these DE-ERGs. The significant enrichment in BP
   (e.g., response to OS, peptide hormone signalling, and regulation of
   defence responses) aligns with pathological features of SONFH,
   including OS and inflammatory responses. GC-induced OS promotes
   apoptosis of osteocytes, primarily by increasing NOX-mediated
   production of ROS^[134]37. Similarly, the exacerbation of inflammatory
   processes by BMP-2–mediated downregulation of IL-34, contributes to
   osteoclast differentiation and the progression of ONFH^[135]38.
   Moreover, the CC enrichment in membrane microdomains and mitochondrial
   outer membranes suggests changes in cellular architecture that could
   influence signal transduction and metabolic processes crucial for cell
   survival and function. KEGG pathway analysis mapped DE-ERGs to pathways
   implicated in SONFH, such as lipid metabolism and atherosclerosis. ERS
   is important in this context because it can be exacerbated by
   intracellular lipid imbalances^[136]39. Dysfunctions in lipid
   metabolism can lead to the accumulation of misfolded proteins in the
   ER, triggering stress responses and ultimately apoptosis and
   bone-tissue degradation^[137]40. Moreover, the atherosclerotic
   processes identified by the pathway analysis indicate that lipid
   accumulation affects cellular health not only directly but also
   indirectly by impairing blood flow, thereby exacerbating the hypoxia
   that leads to osteonecrosis. In addition, identification of the
   NOD-like receptor signalling and FoxO signalling pathways underscores
   the roles of DE-ERGs in modulating immune responses and
   apoptosis^[138]41. This emphasises their potential as therapeutic
   targets and highlights the complex interplay between metabolic
   dysfunction and immune-regulatory mechanisms in SONFH.

   The PPI network analysis provided insight into the molecular
   interactions and pathways linked to the pathophysiology of SONFH. The
   critical clusters in the PPI network, which were related to
   inflammation, immune response, OS, and lipid metabolism, indicates
   steroids affect the activities and interactions of these pathways in
   SONFH. The cytoHubba plugin identified CXCL8, STAT3, IL1B, TLR4, PTGS2,
   TLR2, CASP1, CYBB, CAT, and HMOX1 as hub genes. These genes regulate
   cytokine signalling, immune responses, and stress mechanisms, which are
   important for maintaining cellular homeostasis and responding to
   external stressors including steroids. Their central positions in the
   network implicate them in both normal cellular functions and in
   pathological conditions, potentially driving the processes that lead to
   bone degradation in SONFH. GSEA of hub genes highlighted their effects
   on metabolic regulation, immune modulation, and cellular integrity in
   SONFH. These genes are multifunctional: CXCL8 affects metabolic and
   stress management processes, STAT3 influences cellular growth and
   response to infection, and IL1B integrates metabolic and inflammatory
   pathways, which are key in the pathogenesis of SONFH. TLR4 and PTGS2
   modulate immune responses and apoptosis, suggesting that they have
   potential as therapeutic targets to reduce osteocyte damage. The
   involvement of CASP1 in metabolic processing and cellular clean-up,
   together with the functions of CYBB and TLR2 in immune homeostasis,
   highlight the complex interplay of responses to steroid stress. The
   functions of CAT and HMOX1 in metabolic pathways suggest strategies to
   enhance cellular resilience against SONFH.

   The regulation by miRNAs of the ERS-induced unfolded protein response
   suggests that their effects on pathways that govern cell survival or
   death influence cellular-stress outcomes^[139]42. Indeed, the fact that
   the mRNA-miRNA network centred on hub genes intersecting with ERS
   revealed the importance of miRNAs in SONFH. Among them, miR-23a
   modulates pathways critical for bone health. miR-23a inhibits
   osteogenesis by targeting LRP5, which is crucial for bone homeostasis,
   and icariin stimulates bone formation by modulating this pathway,
   suggesting therapeutic potential^[140]34,[141]35. In addition, miR-23a
   modulates apoptosis and differentiation via LINC00473 in a PLGA
   hydrogel model, suggesting its importance for bone health^[142]36.
   Ganesan et al. found that downregulation of miR-23a led to upregulation
   of TLR2, enhancing the activity of autophagy-related pathways^[143]43.
   Lingam et al. reported that CXCL8 and miR-23a respond to inflammation
   and hypoxia^[144]44. Our predictions suggest that targeting miR-23a
   could improve SONFH outcomes by modulating TLR2 and CXCL8. These hub
   genes, which are significantly modulated by miRNAs, provide insight
   into the complex molecular mechanisms of SONFH and are potential
   targets for miRNA-based therapies.

   Next, we identified potential therapeutic agents that could modulate
   these hub genes. Several drugs, including dibenziodolium and
   N-acetyl-L-cysteine, which have anti-inflammatory and antioxidant
   activity, respectively, potentially interacted with key hub genes.
   Dibenziodolium reduces inflammation-driven exacerbation, and
   N-acetyl-L-cysteine alleviates OS and inflammation, which are important
   in the pathogenesis of SONFH^[145]37. Similarly, simvastatin shows
   therapeutic promise for SONFH because of its anti-inflammatory activity
   and inhibition of the GC receptor^[146]45. Curcumin inhibits M1-type
   macrophage infiltration and osteocyte apoptosis in the femoral head,
   which are linked to the development of SONFH^[147]46. In combination,
   these drugs show greater therapeutic efficacy for SONFH, and may
   improve patient outcomes by modulating the expression of hub genes.

   The significant correlations between immune cells and hub genes suggest
   immune modulation to be a pivotal component of therapeutic strategies.
   The positive correlation between activated mast cells, Tregs, and
   resting dendritic cells suggests that a protective regulatory mechanism
   is disrupted in SONFH. Therapeutic strategies that enhance Treg
   function or stabilise mast-cell activity could mitigate inflammatory
   processes in SONFH. Enhancing the function of Tregs to suppress
   inflammatory responses and pro-inflammatory cells (such as Th17) can
   restore the immune balance and reduce tissue damage, thereby promoting
   the repair and regeneration of bone tissue^[148]47. Conversely, the
   reduction in dendritic cells suggests impairment of immune priming,
   indicating that therapies boosting dendritic cell function could
   restore effective immune responses. He et al. showed that resting
   dendritic cells reduce osteonecrosis inflammation and promote bone
   repair by modulating antigen presentation and cytokine signalling,
   particularly by inhibiting the IL-17 pathway^[149]48.

   Macrophage polarisation into the pro-inflammatory M1 or
   anti-inflammatory M2 phenotype in response to environmental cues is
   pivotal for bone homeostasis, and its dysregulation can lead to
   inflammatory responses that disrupt bone metabolism^[150]49.
   Furthermore, activated NK cells and CD8^+ T cells are components of the
   cell-mediated cytotoxic immune response, and are important for clearing
   damaged cells. However, their overactivation may lead to detrimental
   effects on bone tissue, indicating the need for a balanced immune
   response to maintain skeletal integrity and prevent exacerbation of
   osteonecrosis^[151]50,[152]51. The positive correlations of TLR4, CYBB,
   and CASP1 with M1 macrophages suggest a pro-inflammatory phenotype that
   could exacerbate bone loss and turnover. Targeting these genes to
   modulate M1 macrophage activity may attenuate this pro-inflammatory
   state and induce a shift towards bone preservation and repair.
   Conversely, the negative correlations of TLR4, TLR2, IL1B, HMOX1, CYBB,
   and CASP1 with activated NK cells suggest a protective mechanism
   against the cytotoxicity typically induced by these immune cells.
   Similarly, the negative correlations of these hub genes, together with
   STAT3 and PTGS2, with CD8 T cells suggest a dampening of cytotoxic
   T-cell activity, which could reduce bone tissue damage in SONFH.
   Therefore, adjusting hub gene activity could correct the immune
   imbalance in SONFH, reducing pro-inflammatory effects and fortifying
   defences against overactive immune responses, a promising direction for
   targeted therapy development.

   Finally, we validated the expression levels of 10 hub genes using the
   [153]GSE74089 dataset, confirming the accuracy of our findings from
   [154]GSE123568. To further substantiate these findings and directly
   link them to our study’s focus on ERS-related genes in SONFH, we
   conducted qRT-PCR analysis on peripheral blood samples collected from
   three patients with SONFH and three healthy adults following steroid
   administration. This experimental design aimed to measure the
   expression levels of these hub genes in a clinical context, providing a
   more accurate reflection of their roles in SONFH and their association
   with ERS.

   The qRT-PCR process involved several critical steps: isolating total
   RNA from the blood samples, converting it into cDNA using reverse
   transcription, and then amplifying specific gene targets using
   quantitative PCR to measure their expression levels. This rigorous
   approach ensured that our results were both reliable and relevant to
   the pathophysiological conditions of SONFH, particularly concerning
   ERS. The qRT-PCR analysis yielded consistent results for four hub
   genes: STAT3, IL1B, TLR2, and CASP1. These genes were identified as
   significantly altered in SONFH patients compared to healthy controls,
   underscoring their potential as biomarkers and therapeutic targets.
   Researches have reported that, STAT3 exacerbates SONFH by promoting
   osteoclast differentiation^[155]52, IL1B polymorphisms affect genetic
   susceptibility^[156]53, TLR2 enhances BMSC-mediated angiogenesis and
   osteogenesis^[157]54, and CASP1 affects osteogenic functions via the
   ROS/NLRP3 pathway^[158]55. In an ROC curve analysis, the 10 genes had
   AUC values > 0.6, underscoring their diagnostic potential. This
   indicates that the other hub genes might still be involved in the
   pathophysiology of SONFH, even though their qRT-PCR results were not
   consistent.

   This validation step, by bridging bioinformatics predictions with
   clinical sample analysis, strengthens our findings and confirms the
   relevance of these hub genes in the pathophysiology of SONFH,
   specifically in the context of ERS. Thus, our qRT-PCR results not only
   support the bioinformatic data but also provide a tangible link to
   clinical outcomes and the role of ERS in SONFH, making them crucial for
   developing targeted therapies for this condition.

Strengths and limitations

   This study provides valuable insights into the molecular framework of
   SONFH by identifying and characterizing DE-ERGs. Our comprehensive
   analysis of the [159]GSE123568 dataset revealed significant gene
   expression changes associated with steroid exposure, offering a robust
   foundation for understanding the cellular mechanisms involved in SONFH.
   The use of GO and KEGG enrichment analyses, PPI network analysis, and
   validation with the [160]GSE74089 dataset further strengthens the
   reliability of our findings.

   However, the study has limitations. Firstly, the lack of randomization
   may introduce selection bias and affect the internal validity of our
   findings. Future studies should include randomized controlled trials to
   minimize bias and enhance the robustness of the results. Additionally,
   the generalizability of our findings is constrained by the specific
   characteristics of our study population, which was drawn from a single
   institution. This may not fully represent the broader population of
   patients with SONFH. Future studies should include more diverse and
   larger sample sizes from multiple centers to improve the applicability
   of the findings. The small sample size for qRT-PCR analysis may also
   limit the generalizability of the results. Future research should
   include in vitro validation of the predicted miRNA regulatory
   mechanisms and the therapeutic potential of the identified drugs to
   gain deeper insights into their clinical implications.

Conclusion

   Using the [161]GSE123568 dataset, we identified 195 DE-ERGs in SONFH,
   which were related to oxidative stress, immune responses, and metabolic
   pathways. GO and KEGG analyses, and the PPI network, identified CXCL8,
   STAT3, IL1B, TLR4, PTGS2, TLR2, CASP1, CYBB, CAT, and HMOX1 as key hub
   genes, which regulate metabolic, immune, and cellular integrity
   processes. The identification of 261 miRNAs, particularly hsa-miR-23,
   emphasised the complexity of the regulation of genes that modulate bone
   metabolism. Moreover, we found that dibenziodolium,
   N-acetyl-L-cysteine, simvastatin, and curcumin target inflammatory
   pathways important in the pathophysiology of SONFH. Twenty-two immune
   cell subtypes showed significant interactions, notably activated mast
   cells and Tregs, and there was a reduction in the population of
   dendritic cells, indicating altered immune function in SONFH. In
   addition, the positive and negative correlations of hub genes with
   immune cell types underscore their potential as therapeutic targets.
   Validation by cross-dataset comparison, qRT-PCR and ROC analysis
   confirmed the upregulation of STAT3, IL1B, TLR2, and CASP1. Further
   validation by ROC curve analysis confirmed that these hub genes are
   reliable diagnostic markers. The findings provide molecular insight
   into SONFH and will facilitate the development of targeted
   interventions.

Acknowledgements