Identification of Ferroptosis-Related Biomarkers for Diagnosis and Molecular Classification of Staphylococcus aureus-Induced Osteomyelitis

Abstract

Objective

Staphylococcus aureus (SA)-induced osteomyelitis (OM) is one of the most common refractory diseases in orthopedics. Early diagnosis is beneficial to improve the prognosis of patients. Ferroptosis plays a key role in inflammation and immune response, while the mechanism of ferroptosis-related genes (FRGs) in SA-induced OM is still unclear. The purpose of this study was to determine the role of ferroptosis-related genes in the diagnosis, molecular classification and immune infiltration of SA-induced OM by bioinformatics.

Methods

Datasets related to SA-induced OM and ferroptosis were collected from the Gene Expression Omnibus (GEO) and ferroptosis databases, respectively. The least absolute shrinkage and selection operator (LASSO) and support vector machine-recursive feature elimination (SVM-RFE) algorithms were combined to screen out differentially expressed-FRGs (DE-FRGs) with diagnostic characteristics, and gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) were used to explore specific biological functions and pathways. Based on these key DE-FRGs, a diagnostic model was established, and molecular subtypes were divided to explore the changes in the immune microenvironment between molecular subtypes.

Results

A total of 41 DE-FRGs were identified. After screening and intersecting with LASSO and SVM-RFE algorithms, 8 key DE-FRGs with diagnostic characteristics were obtained, which may regulate the pathogenesis of OM through the immune response and amino acid metabolism. The ROC curve indicated that the 8 DE-FRGs had excellent diagnostic ability for SA-induced OM (AUC=0.993). Two different molecular subtypes (subtype 1 and subtype 2) were identified by unsupervised cluster analysis. The CIBERSORT analysis revealed that the subtype 1 OM had higher immune cell infiltration rates, mainly in T cells CD4 memory resting, macrophages M0, macrophages M2, dendritic cells resting, and dendritic cells activated.

Conclusion

We developed a diagnostic model related to ferroptosis and molecular subtypes significantly related to immune infiltration, which may provide a novel insight for exploring the pathogenesis and immunotherapy of SA-induced OM.

Keywords:

• osteomyelitis
• ferroptosis
• molecular subtype
• immune infiltration
• biomarker

Introduction

Osteomyelitis (OM) is a complex bone disease that causes inflammation and the destruction of osseous tissue. It usually occurs in the long bones, especially femur, tibia, and humerus.^1,2 Three clinical mechanisms lead to bone infection: OM resulting from the spread of hematogenous sources, OM occurring secondary to adjacent soft tissues and joints, and direct inoculation of bacteria into the bone as a result of trauma or surgery.³ The Gram-positive bacterium Staphylococcus aureus (SA) is the most implicated microorganism.⁴ Since the availability of antibiotics, mortality rates from OM, including SA-induced OM, have decreased significantly, but the infection rate and death rate are still significant.^5,6 The key to successful management is early diagnosis, including bone sampling for microbiological and pathological examination to allow targeted and long-lasting antimicrobial therapy. Therefore, screening and searching for key biomarkers are crucial for the early detection of SA-induced OM, which is expected to significantly improve the healing of patients.

Cell death is critical for the development of multicellular organisms, and necrosis and apoptosis are the main forms of cell death.⁷ In recent years, nonapoptotic cell death has been proven to be cell death under gene regulation, while iron death is a unique form of nonapoptotic cell death.⁸ In 2012, Dixon et al⁹ first proposed the concept of ferroptosis, a form of regulated cell death, is defined as an iron-catalyzed form of regulated necrosis that is characterized by mitochondrial atrophy and increased mitochondrial membrane density, the accumulation of iron and lipid reactive oxygen species and the involvement of a unique set of genes, which is under the constitutive control of glutathione peroxidase 4, ferroptosis through excessive peroxidation of polyunsaturated fatty acids.^10–12 In addition, the relationship between ferroptosis and inflammatory diseases has been investigated by some groups recently and they verified the positive anti-inflammatory effect of ferroptosis in inflammation, including nonalcoholic steatohepatitis,^13,14 lung fibrosis,¹⁵ and ischemia-reperfusion injury.^16,17 The pathophysiological process of OM involves immune and inflammatory reactions, and ferroptosis with immunogenicity may play an active role in the pathogenesis and treatment of OM.^18,19 However, the mechanism of ferroptosis in SA-induced OM is still unclear.

In this study, based on the Gene Expression Omnibus (GEO) database, we used machine learning algorithms to screen the optimal differentially expressed ferroptosis-related genes (DE-FRGs) with diagnostic characteristics. Receiver operating characteristic (ROC) curves were performed to evaluate the diagnostic ability of these DE-FRGs, and an independent dataset was validated. Then, the nomogram model was constructed to predict the incidence rate of SA-induced OM. Additionally, based on these key DE-FRGs, further unsupervised cluster analysis was used to divide osteomyelitis samples into different molecular subtypes to explore the changes in the abundance of immune cell infiltration among the different molecular subtypes.

Materials and Methods

Data Acquisition and Preprocessing

In this study, all gene expression datasets (GSE6269, GSE16129 and GSE30119) for SA-induced OM and normal samples were downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database. After eliminating the batch effect²⁰ and excluding the “bacteraemia” samples, the GSE6269 and GSE16129 datasets were combined as the training set for this study, including 85 OM samples and 35 normal samples. The GSE30119 datasets were extracted as the test set, including 39 OM and 44 healthy samples. The conversion of gene symbols and the merging of data were performed using the “perl” software.²¹ Additionally, 348 FRGs were preprepared and obtained from the FerrDb database²² for further analysis (Supplementary Table S1).

Identification of Differentially Expressed Genes (DEGs)

Based on the 348 FRGs, FRGs were extracted from 85 OM and 35 healthy samples in the training set. To compare the DE-FRGs between OM and normal samples, differential expression analysis was performed using the “limma” package of R software²³ and “pheatmap”²⁴ package were employed for visualization. The screening criterion for DE-FRGs was P value<0.05. To further explore the correlation between these DE-FRGs, a correlation diagram was constructed using the “corrplot” package.

Functional Enrichment Analysis

To explore the biological functions of these DE-FRDs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were employed to search for significantly enriched biological functions and pathways using the “clusterProfiler” package.²⁵ GO functional analysis included the top 10 significant biological processes (BP), cellular components (CC) and molecular functions (MF), and KEGG pathway analysis included the top 20 significant pathways. The criteria for screening were set as a P value<0.05 and adjusted P value<1.

Identification of Optimal Diagnostic Biomarkers and Construction of Nomogram Model for SA-Induced OM

To screen the optimal gene biomarkers to distinguish osteomyelitis and normal samples, the least absolute shrinkage and selection operator (LASSO) algorithm was employed with the “glmnet” package to reduce the dimensionality of the data, and to further identify significant gene biomarkers of osteomyelitis based on DE-FRGs. Meanwhile, a support vector machine-recursive feature elimination (SVM-RFE) model was constructed using the “e1071” package to rank the importance of the significant genes, followed by 10-fold cross-validation to obtain the error and accuracy of cross-validation. Additionally, the LASSO and SVM-RFE algorithms were combined to identify the optimal gene biomarkers for osteomyelitis, and the results of the cross-validation were visualized using the “VennDiagram” package. Receiver operating characteristic (ROC) curves were constructed to assess the ability of the optimal gene biomarkers to differentiate between osteomyelitis and normal samples, and the area under the curve (AUC) and accuracy were measured. Finally, logistic regression models were constructed using the ‘glmnet’ package to predict whether samples in the training set belonged to osteomyelitis samples or normal samples, and a ROC curve was plotted to assess the overall diagnostic ability.

Furthermore, to explore the specific role of key DE-FRGs in the diagnosis of OM, the “rms” package was used to construct a nomogram model to predict the incidence of osteomyelitis. Calibration curves and decision curve analysis (DCA) were used in the analysis to assess the fit of the model.²⁶

Single-Gene Gene Set Enrichment Analysis (GSEA) Enrichment Analysis

To further analyze the related pathways and potential biological functions of these optimal DE-FRGs in OM, we used the “enrichplot” and “clusterProfiler” packages to perform GSEA enrichment analysis for each signature gene, with two gene sets “c5.go.symbols.gmt” and “c2.cp.kegg.symbols.gmt” as the predefined sets. The top 6 pathways with significant enrichment were visualized and the screening threshold was set at P value<0.05.

Single-Gene Gene Set Variation Analysis (GSVA) Enrichment Analysis

GSVA is a method to estimate the variation of pathway activity in samples with an unsupervised way due to its stability and is often used in the data analysis of gene expression profiles.²⁷ Each signature gene was analyzed by GSVA based on predefined sets “c5.go.symbols.gmt” and ‘c2.cp.kegg.symbols.gmt’. First, the samples are scored and corrected. Then, the samples were divided into groups according to the expression of the target gene, and the difference in GSVA score between the high expression group and the low expression group was further analyzed. The screening condition for significant differences was |t| >2 and P value<0.05.

Immune Microenvironment of SA-Induced OM

The Cell type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) is a method to characterize the cellular composition of tissues through gene expression profiling and is used to find cellular biomarkers and therapeutic targets.²⁸ CIBERSORT algorithm was used to estimate the proportion of 22 immune cell types in the training set. Subsequently, differences in the abundance of 22 immune cell infiltrates between normal and osteomyelitis samples were observed and a violin plot was constructed using the “vioplot” package.²⁹ Finally, to explore the relationship between the signature genes and the immune microenvironment, a correlation analysis was performed between all the signature genes and the 22 immune cells.

Identification of Molecular Subtypes Based on Key DE-FRGs

In consensus clustering, the K-means algorithm is run multiple times to obtain the input partition, and the common matrix is calculated based on the partition result. The ultimate goal is to detect sample subtypes with similar characteristics.³⁰ The “ConsensusClusterPlus” package of R software was used to cluster osteomyelitis samples with significant DE-FRGs to create subtypes with differential characteristics.^25,31 The most appropriate number of clusters was screened based on the consistency heatmap, consistency scores, cumulative distribution function (CDF) and CDF delta area curve. Specifically, the value of K corresponding to when the CDF reaches its approximate maximum value is the best grouping result.³² The maximum value of the number of clusters K was set to 10. Principal component analysis (PCA) was performed on the classified samples to determine whether the groupings accurately reflected the characteristics of OM patients.³³ Furthermore, the expression levels of significant DE-FRGs were further analyzed to determine whether they differed significantly between molecular subtypes.

Verification of Diagnosis and Expression Based on the Validation Set

To verify the distinguishing ability of key DE-FRGs between OM samples and healthy samples, the expression levels of DE-FRGs in the GSE30119 dataset were extracted by the “limma” package, and the Logistic regression model and corresponding ROC curve were constructed to evaluate the overall diagnostic ability. In addition, based on the expression of significant DE-FRGs, we validated whether the expression level of each significant DE-FRG significantly differed between 39 SA-induced OM samples and 44 healthy samples.

Clinical Correlation Analysis of the Validation Set

To further explore the clinical application value of the significant DE-FRGs, the age characteristics of 39 SA-induced OM samples and 44 healthy samples in the validation set were extracted. The age correlation of the significant DE-FRGs was analyzed one by one, and the clinical correlation diagrams were drawn. A correlation coefficient greater than 0 indicates a positive correlation and a correlation coefficient less than 0 indicates a negative correlation.

Establishment of SA-Induced OM Model in Rats and qPCR Validation of the Bioinformatics Results

36 Sprague‒Dawley rats weighing 500 g were purchased from Kunming Chushang Co., Ltd. (Kunming, China), of which 18 were in the experimental group and 18 in the control group. The experimental group was injected with 6 μL of 1×10⁶ SA at the tibia site, and the control group was injected with sterile saline at the same site. Sodium pentobarbital was injected intraperitoneally for anesthesia. For the specific procedure of OM model establishment, please refer to Zak et al³⁴ and our previously published article.³⁵ All experimental procedures were approved by the 920th Hospital of Joint Logistics Support Force Committee on Ethics for the care and use of laboratory animals (2022–075-01) and conducted in accordance with the Guide for the Care and Use of Laboratory Animals (National Institutes of Health, USA).

Four weeks after the establishment of the OM model, total RNA was extracted from the tibial lesion tissues of six rats in the experimental group and from the tibial tissues of six rats in the control group. Reverse transcription of mRNA using SweScript RT I First strand cDNA synthesis kit (Service Bio, Guangzhou, China). After a short centrifugation, the reaction was performed for 40 cycles at the CFX96 real-time quantitative fluorescence PCR instrument (Bio-Rad, Hercules, CA, USA) under the following conditions: pre-denaturation 95 °C, 1 min; denaturation 95 °C, 20s; annealing 55 °C, 20s; and extension 72 °C, 30s. The ethnographic primer sets used were as follows: SLC38A1: forward 5′-TCCTTCAGCCATAAAATCCCT-3′ and reverse 5′-TACCATCACCACCAACACTCG-3′; MAPK9: forward 5′-GAATCCGAACGAGACAAAAT-3′ and reverse 5′-CGGGGTCATACCAAACAGTA-3′; SNCA: forward 5′-GCGTCCTCTATGTAGGTTCC-3′ and reverse 5′-GTTCTTTGGTCTTCTCAGCC-3′; KLF2: forward 5′-CTATCTTGCCGTCCTTTGCC-3′ and reverse 5′-TGTTTAGGTCCTCATCCGTG-3′; EGR1: forward 5′-TAAAGGCAATGACAAGGGAG-3′ and reverse 5′-TCAGGCAAGTGAGCGTAGAA-3′; STAT3: forward 5′-CCCCACTGGTCTACCTCTAC-3′ and reverse 5′-CAATACTTTCCGAATGCCTC-3′; SREBF2: forward 5′-TGGCAAATCAGAAAAACAAGC-3′ and reverse 5′-TCAGAGTCAATGGAATAGGGG-3′; ABCC5: forward 5′-TGGACGAGGAACATACGAA-3′ and reverse 5′-AGGCACACGATGGACAGGA-3′; GAPDH: forward 5′-CCCATCACCATCTTCCAGG-3′ and reverse 5′-CATCACGCCACAGTTTCCC-3′.

Statistical Analysis

Student’s t-test was used for comparisons between the two groups. Pearson’s test was used to reveal the relationship between the DE-FRGs and immune infiltration. P<0.05 was considered significant. All data analyses were performed in R software (version 4.1.3).

Results

Identification of DE-FRGs in the Training Set

In 85 SA-induced OM samples and 35 normal samples after combining and eliminating batch effects, 41 DE-FRGs were obtained, of which 30 were significantly upregulated and 11 were significantly downregulated (Supplementary Table S2). The results of the clustering heatmap (Figure 1a) showed the expression of 41 DE-FRGs in different samples. Correlation analysis of these genes (Figure 1b) revealed that SLC38A1 was significantly negatively correlated with ATG7, TLR4, CTSB and LAMP2, and significantly positively correlated with CIRBP. MAPK9 showed no correlation with the remaining 40 DE-FRGs, except for one negative correlation with BID.

Figure 1 Expression levels and correlations of DE-FRGs in OM. (a) Differential heatmap of these DE-FRGs. Red indicates the DE-FRG is highly expressed in the sample, and blue indicates the DE-FRG is lowly expressed in the sample. (b) The correlation analysis of these DE-FRGs. Red represents positive correlation, blue represents negative correlation. *: P<0.05, **: P<0.01, ***: P<0.001.

Functional Enrichment Analysis of DE-FRGs

We further explored the biological functions and pathways of 41 DE-FRGs by gene enrichment analysis. In the GO enrichment analysis (Figures 2a and b), response to oxidative stress, ficolin−1−rich granule and DNA−binding transcription activator activity, and RNA polymerase II−specific were the most significantly enriched GO terms for BP, CC, and BF, respectively. In the KEGG pathway analysis (Figure 2c), DEGs were mainly involved in necroptosis, lipid and atherosclerosis, and HIF−1 signaling pathway pathways. Consequently, the above evidence suggested that DE-FRGs may play a role in the pathogenesis of SA-induced OM by participating in the regulation of apoptosis and transcription factors.

Figure 2 Functional enrichment analysis of GO and KEGG pathway in DE-FRGs. (a) GO enrichment analysis of top 20 biological processes. (b) Top 10 enriched GO terms for biological processes (BP), cellular components (CC), and biological functions (BF). (c) Enrichment analysis of the top 20 KEGG pathways.

Diagnostic Biomarkers Were Identified for OM Based on Significant DE-FRGs

To distinguish between patients with OM and healthy samples, we further screened the diagnostic genes, and introduced two machine learning algorithms, LASSO and SVM-RFE, into the training sets (GSE6269 and GSE16129 datasets). Specifically, Logistic regression algorithm was used for 10 rounds of cross-validation, and 11 significant gene biomarkers were selected from 41 DE-FRGs (Figure 3a and b). Then, 41 DE-FRGs were filtered by the SVM-RFE algorithm, and 18 significant gene biomarkers were selected (Figure 3c and d). The genetic biomarkers obtained from the LASSO and SVM-RFE models were intersected, and 8 significant biomarkers were finally obtained (SLC38A1, MAPK9, SNCA, KLF2, EGR1, STAT3, SREBF2 and ABCC5) (Figure 3e). The diagnostic curve of ROC was drawn based on 8 significant biomarkers, and the AUC values of all genes were greater than 0.65, indicating that each significant gene biomarker can well distinguish OM samples from healthy samples (Figure 4a). Interestingly, the “glmnet” package was used to establish a Logistic regression model, and the results of the ROC curve showed that the regression model based on eight characteristic gene markers had higher accuracy than the diagnostic efficiency of a single gene marker (AUC= 0.993, 95% CI: 0.979–0.100) (Figure 4b and Supplementary Table S3). Finally, the nomogram model was constructed by using the “rms” package to predict the incidence of OM (Figure 4c). The results of the calibration curve revealed that the prediction of the nomogram model was basically consistent with that of ideal model (Figure 4d), and the clinical decision based on this model could be beneficial to patients with SA-induced OM (Figure 4e).

Figure 3 8 Diagnostic biomarkers were identified for OM based on significant DE-FRGs. (a) Coefficient distribution graph; a line represents a gene, and the ordinate of each gene corresponds to a coefficient. (b) Screening of gene biomarkers from DE-FRGs using the LASSO algorithm. (c) Result of cross-validation error. (d) Result of cross-validation accuracy. (e) Gene biomarkers are obtained based on the intersection of LASSO and SVM-RFE algorithms.

Figure 4 ROC curve and nomogram model based on 8 diagnostic biomarkers. (a) ROC curves for the 8 diagnostic biomarkers. (b) a Logistic regression model was constructed to identify the OM samples. (c) Nomogram of 8 gene biomarkers in the diagnosis of OM patients. (d) Calibration curve. (e) Decision curve analysis of the nomogram model.

Diagnostic Biomarkers Were Associated with OM-Related Enrichment Pathways

To further explore the molecular mechanism of 8 characteristic biomarkers related to the diagnosis of SA-induced OM, each gene biomarker was analyzed by ssGSEA-KEGG pathway enrichment analysis. The figures showed the top 6 enrichment pathways at most (Figures 5a-d and Supplementary Figure S1). The results of comprehensive analysis indicated that 8 characteristic biomarkers were significantly enriched in the cell cycle, fc-gamma-r mediated phagocytosis, lysosome, ribosome, spliceosome, axon guidance, cell adhesion molecules cams, receptor interaction (ECM receptor interaction, neuroactive ligand receptor interaction and cytokine-cytokine receptor interaction) and lipid metabolism (glycerophospholipid metabolism and inositol phosphate metabolism). In addition, these gene biomarkers were also significantly enriched in multiple signaling pathways (nod like receptor signaling pathway, TGF-beta signaling pathway and ERBB signaling pathway) and a variety of diseases (small cell lung cancer, arrhythmogenic right ventricular cardiomyopathy and amyotrophic lateral sclerosis).

Figure 5 Single-gene GSEA-KEGG pathway analysis in SLC38A1 (a), MAPK9 (b), SNCA (c), KLF2 (d).

Based on the median expression level of each biomarker, SA-induced OM samples were divided into a high-expression group and a low-expression group, and the enrichment pathways between the two groups were explored by GSVA enrichment analysis. The results of comprehensive analysis suggested that high expression of SLC38A1 may activate complement and coagulation cascades, pantothenate, and CoA biosynthesis and citrate cycle TCA cycle pathways, while low expression of SLC38A1 may activate nicotinate and nicotinamide metabolism and proximal tubule bicarbonate reclamation pathways (Figure 6a). High expression of EGR1 was associated with alpha-linolenic acid metabolism, DNA replication, and arachidonic acid metabolism, while low expression of EGR1 was associated with the intestinal immune network for IgA production and glycosphingolipid biosynthesis lacto and neolacto series (Figure 6b). Consistent with EGR1, high expression of KLF2 may also activate arachidonic acid metabolism and alpha-linolenic acid metabolism pathways, while low expression of KLF2 may activate autoimmune thyroid disease and proximal tubule bicarbonate reclamation pathways. In the STAT3 high-expression group, amino acid metabolism-associated pathways (tyrosine metabolism, nicotinate and nicotinamide metabolism, and butanoate metabolism) were significantly enriched, while complement and coagulation cascades and glyoxylate and dicarboxylate metabolism were enriched in the STAT3 low-expression group. High expression of SREBF2 was only related to retinol metabolism and linoleic acid metabolism, while low expression of SREBF2 was mainly associated with the intestinal immune network for IgA production and glycosphingolipid biosynthesis lacto and neolacto series. Similarly, high expression of SNCA was only related to nicotinate and nicotinamide metabolism and butanoate metabolism, while low expression of SNCA was mainly associated with glyoxylate and dicarboxylate metabolism and glycine serine and threonine metabolism. The highly expressed MAPK9 was only related to the base excision repair pathway and highly expressed ABCC5 was only related to taurine and hypotaurine metabolism (Supplementary Figures S2). These results suggested that 8 significant biomarkers may regulate the pathogenesis of SA-induced OM through immune response- and amino acid metabolic-associated pathways.

Figure 6 Differentially activated pathways between the high- and low-expression groups based on the expression levels of each gene biomarkers. (a) GSVA-KEGG pathway analysis in SLC38A1. (b) GSVA-KEGG pathway analysis in EGR1.

Immune Microenvironment Analysis

The pathophysiological mechanism of SA-induced OM involves inflammation and immune response, and many immune cells and factors are involved in its pathogenesis.³⁶ Therefore, to explore the immune response mechanism of SA-induced OM, we used the CIBERSORT algorithm to evaluate the abundant difference of immune cells between OM patients and healthy patients (Figure 7a). The results showed that the abundance of T cells follicular helper, T cells regulatory (Tregs), monocytes, macrophages M0, macrophages M2, mast cells resting, eosinophils and neutrophils in OM samples was significantly higher than that in control samples, while B cells memory, T cells CD4 naive and NK cells resting were significantly higher than those in control samples. The results of immune correlation analysis indicated that both SLC38A1 and KLF2 had a strong negative correlation with macrophages M2, while SREBF2, STAT3 and EGR1 had a strong positive correlation with neutrophils (Figure 7b). The immune microenvironment of SA-induced OM may be related to SLC38A1, KLF2, SREBF2, STAT3 and EGR1.

Figure 7 Immune microenvironment analysis. (a) Implemented the CIBERSORT algorithm to evaluate the abundant difference of immune cells between OM samples and healthy samples. (b) The results of immune correlation analysis indicated that both SLC38A1 and KLF2 had a strong negative correlation with macrophages M2, while SREBF2, STAT3 and EGR1 had a strong positive correlation with neutrophils. *P<0.05, **P<0.01, ***P<0.001.

Results of Molecular Subtypes Based on Significant DE-FRGs

Based on 8 significant biomarkers, a total of 9 cluster results were established. To screen the most suitable molecular subtype, the consistency score and CDF value were as used references (Figures 8a-c). The optimal clustering result was K=2, and the consistency score was close to 1.0. To verify the differences between subtype 1 and subtype 2 samples, PCA showed that there were significant differences in transcriptome among the different subtypes (Figure 8d). Then, the results of the boxplot and correlation heatmap showed that MAPK9, SNCA, KLF2 and STAT3 were highly expressed in subtype 1, while SLC38A, EGR1 and ABCC5 were highly expressed in subtype 2 (Figures 9a and b). Finally, we explored the differences in immune microenvironment characteristics between subtype 1 and subtype 2 (Figures 9c and d). The results showed that subtype 1 OM had higher immune cell infiltration rates, mainly in T cells CD4 memory resting, macrophages M0, macrophages M2, dendritic cells resting, and dendritic cells activated. T cells CD4 naive and NK cells resting had higher infiltration abundance in subtype 2 OM.

Figure 8 Molecular subgroups based on clustering analysis of 8 gene biomarkers. (a) Heatmap of 2 clusters (k = 2) based on gene biomarkers. (b) Consistency scores of 2–9 clusters. (c) Cumulative distribution graph. (d) PCA analysis of the 2 clusters: blue indicates cluster 1 samples; red indicates cluster 2 samples.

Figure 9 The differences in gene expression levels and immune microenvironment characteristics between two different clusters. (a) Differences in the expression levels of 8 gene biomarkers between cluster 1 and cluster 2: blue indicates cluster 1 samples; red indicates cluster 2 samples. (b) Clusting heatmap of gene expression levels between two different clusters: red represents high expression, blue represents low expression. (c) Immune cell content stacking plot between cluster 1 and cluster 2: different colors indicate different immune cells; the horizontal axis is the different clusters. (d) Immune cell content histogram: the horizontal axis indicates 22 immune cells; the vertical axis indicates infiltration abundance; blue indicates cluster 1 samples; red indicates cluster 2 samples. *: P<0.05, **: P<0.01, ***: P<0.001.

Validation of Diagnosis, Expression and Clinical Correlation Analysis Based on the Validation Set

To verify the ability of the Logistic regression model based on 8 significant biomarkers to distinguish SA-induced OM from normal samples, the GSE30119 dataset was used as the validation set. The results of the ROC curve analysis suggest that 8 characteristic significant biomarkers distinguish OM from normal samples (Figure 10a). In addition, we verified the difference in the expression of 8 characteristic biomarkers between 39 OM samples and 44 healthy samples (Figures 10b-i), suggesting that SNCA, EGR1, SREBF2, MAPK9, and STAT3 were highly expressed and SLC38A1, KLF2, and ABCC5 expressed at low levels in OM samples (P<0.05), which is consistent with the results of the training set analysis. The results of clinical correlation analysis showed that except for the positive correlation between STAT3 and the ages of patients with OM, the other gene biomarkers had no significant correlation with age (Figures 11a-h).

Figure 10 Logistic regression model and expression levels of gene biomakers in the validation set. (a) Logistic regression model in the validation set. Expression levels of SLC38A1 (b), MAPK9 (c), SNCA (d), KLF2 (e), EGR1 (f), STAT3 (g), SREBF2 (h) and ABCC5 (i) in the validation set.

Figure 11 Clinical correlation analysis between 8 gene biomakers and the ages of OM patients. Clinical correlation analysis of SLC38A1 (a), MAPK9 (b), SNCA (c), KLF2 (d), EGR1 (e), STAT3 (f), SREBF2 (g) and ABCC5 (h) in the validation set: R>0 indicates the two variables are positively correlated; R<0 indicates the two variables are negatively correlated; P<0.05 represents that the two variables are significantly correlated.

RT-qPCR Validation of the Expression Levels of 8 Gene Biomarkers

To verify the expression levels of the 8 gene markers, rat model of OM was established and validated by RT-qPCR experiments (Figures 12a-h). RT-qPCR results suggested that MAPK9, SNCA, EGR1, STAT3, and SREBF2 were highly expressed in OM lesion tissues, while SLC38A1 and KLF2 were highly expressed in the control group. The trends of gene expression were fully consistent with those of the training set (GSE6269 and GSE16129 datasets) and the validation set (GSE30119 dataset). Additionally, although there was no significant difference in mRNA expression level of ABCC5 between OM group and control group, it was consistent with the trend of gene expression in GSE30119 dataset.

Figure 12 RT-qPCR validation. Rat models of OM was established and total RNA was extracted from the focal bone tissue and healthy bone tissue of rat model of OM, the mRNA expression levels of SLC38A1 (a), MAPK9 (b), SNCA (c), KLF2 (d), EGR1 (e), STAT3 (f), SREBF2 (g) and ABCC5 (h) in the the focal bone tissue and control group. ns: P≥0.05, *: P<0.05, **: P<0.01, ***: P<0.001.

Discussion

OM is a persistent inflammatory bone disease, and its diagnosis requires a multidisciplinary approach including recognition and assessment of clinical symptoms, imaging tests and laboratory investigations.^37–39 The gold standard for the diagnosis of SA-induced OM is bone biopsy, which is invasive and expensive, and is often performed after clinical symptoms.⁴⁰ Delayed diagnosis is prone to various serious complications, including malformation and sepsis, and it is difficult to cure in the clinic.⁴¹ In this study, we identified DE-FRGs between SA-induced OM samples and healthy samples by differential expression analysis combined with machine learning to screen the most appropriate features of DE-FRGs. Then, ssGSEA and ssGSVA analyses were used to explore the biological functions and pathways related to characteristic DE-FRGs. In addition, diagnostic models and molecular subtypes of OM were constructed based on 8 characteristic DE-FRGs to further explore the correlation between these DE-FRGs and changes in the immune microenvironment. The results of the training set and validation set proved that the diagnostic model of OM based on 8 key gene biomarkers had good diagnostic efficacy.

The diagnostic model included 8 significant biomarkers (SLC38A1, MAPK9, SNCA, KLF2, EGR1, STAT3, SREBF2 and ABCC5), and the AUC values represented by the area under the ROC curves of all genes were greater than 0.65, indicating that the 8 significant biomarkers have certain accuracy in distinguishing OM samples from healthy samples. Among them, the AUC values of SLC38A1, SNCA, KLF2, EGR1 and SREBF2 were all greater than 0.80. Combined with the results of the correlation with changes in the immune microenvironment, SLC38A1, KLF2, EGR1 and SREBF2 may be closely related to the diagnosis and immunotherapy of SA-induced OM. As a member of the Solute Carrier Family 38 (SLC38), SLC38A1 is mainly expressed in the brain and placenta, and is responsible for the transport of short chain neutral amino acids.^42,43 Numerous studies have confirmed that overexpression of SLC38A1 is associated with poor prognosis, proliferation and migration of many tumors, such as colorectal cancer,⁴⁴ liver cancer,⁴⁵ gastric cancer⁴⁶ and osteosarcoma.⁴⁷ Additionally, SLC38A1 may participate in the regulation of musculoskeletal diseases. Giresi et al⁴⁸ found that the expression level of SLC38A1 was significantly downregulated in patients with sarcopenia, which is similar to our results. In the differential expression analysis results of the training set and test set, SLC38A1 was significantly expressed at low levels with OM. KLF2 gene is mainly involved in the regulation of cell growth and differentiation and has been proven to be able to target Runt related transcription factor 2 (RUNX2) to regulate osteoblast differentiation.⁴⁹ On the other hand, Dipranjan et al⁵⁰ found that KLF2 gene can also regulate the formation of osteoclasts through autophagy. In arthritis model, KLF2 can promote osteoclast differentiation and increase proinflammatory factors.⁵¹ Interestingly, some studies have proven that KLF2 participates in the regulation of macrophages in muscle injury, which is manifested in the increase in the number of macrophages and the enhancement of phagocytosis.⁵² In this study, KLF2 and macrophages M2 showed a significant negative correlation in the immune environment of OM.

EGR1 gene plays an important role in regulating the inflammatory response and is the first transcription factor to fight against inflammatory enhancers.⁵³ A study showed that EGR1 is significantly increased in the cartilage of patients with osteoarthritis, further revealing that EGR1 promotes cartilage deformation through the β-catenin signaling pathway.⁵⁴ Similarly, this study also found that EGR1 was significantly increased in OM samples, indicating that EGR1 may play an important role in the regulatory mechanism of inflammatory diseases. SREBF2 is a transcription factor that regulates lipid homeostasis and participates in the regulation of inflammatory diseases.⁵⁵ A study has confirmed that the upregulation of SREBF2 is involved in the development of osteoarthritis,⁵⁶ which is consistent with our findings.

Even so, our study still has some limitations. First, although GSE30119 was included in this study as a verification set of diagnostic ability, a larger sample size dataset and experimental verification are still needed. Second, this study was only limited to the mRNA level of SA-induced OM samples. Future research should explore the pathological development and early diagnosis of SA-induced OM from multiple histology (metabonomics and lipoomics).

Conclusion

In conclusion, SLC38A1, MAPK9, SNCA, KLF2, EGR1, STAT3, SREBF2 and ABCC5 were identified as characteristic biomarkers. We established a diagnostic model of SA-induced OM patients, which has good diagnostic value. Based on 8 characteristic ferroptosis-related genes, there were significant differences in the immune microenvironment between subtype 1 and subtype 2 OM. This study will contribute to the development of early diagnostic indicators and immunotherapy methods for SA-induced OM.

Abbreviations

OM, Osteomyelitis; SA, Staphylococcus aureus; GEO, Gene Expression Omnibus; LASSO, Least absolute shrinkage and selection operator; SVM-RFE, Support vector machine-recursive feature elimination; GSEA, Gene set enrichment analysis; GSVA, Gene set variation analysis; DEGs, Differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, Biological processes; CC, Cellular components; MF, Molecular functions; ROC, Receiver operating characteristic; AUC, Area under the curve; DCA, Decision curve analysis; CIBERSORT, Cell Type Identification by Estimating Relative Subsets of RNA Transcripts; CDF, Cumulative distribution function; PCA, Principal component analysis; RT-qPCR, Reverse transcription quantitative PCR.

Data Sharing Statement

All the results are presented in the article and Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics Statement

All experimental procedures were approved by the 920th Hospital of Joint Logistics Support Force Committee on Ethics for the care and use of laboratory animals (2022-075-01) and conducted in accordance with the Guide for the Care and Use of Laboratory Animals (National Institutes of Health, USA).

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Disclosure

All the authors declare that they have no competing interests in this work.

Acknowledgments

We are grateful to the authors who gave the GEO public dataset.

Additional information

Funding

This study was funded by National Natural Science Foundation of China (Grant No.81772367, 82072392); the Yunnan Traumatology and Orthopedics Clinical Medical Center (Grant No. ZX20191001); the Grants from Yunnan Orthopedics and Sports Rehabilitation Clinical Medicine Research Center (Grant No. 202102AA310068).