We downloaded the S. aureus infection-related osteomyelitis chip data and the corresponding clinical information for the S. aureus infection-related osteomyelitis dataset GSE30119 from the Gene Expression Omnibus database (Toufiq et al., 2000). The sample source was Homo sapiens, and the sequencing platform used was GPL6947 Illumina HumanHT-12 V3.0. The dataset included two experimental sets (C1 and C2) representing 44 normal samples and 99 samples from patients with S. aureus infections. We included C1 as a training set with 22 normal samples and 40 patients with S. aureus infection and C2 as a validation set with 22 normal samples and 59 patients with S. aureus infection. In order to validate the model, we collected the bone marrow samples of S. aureus infection GSE16129 as the validation set. The R package “sva” was used to correct the batch effect between different datasets and perform log2 standardization (Leek et al., 2012). The expression distributions before and after standardization and batch correction were visualized using a box diagram (Figure 1).

To analyze the effect of m6A-related genes on disease, the R package “limma” was used to perform differential gene analysis on normal and diseased samples in the dataset (Ritchie et al., 2015). The volcano map showed a log2fold change (logFC) absolute value >2 and Padj <0.05, which were set as the differential expression results, and the heatmap showed the top 50 differentially expressed genes between normal samples and disease samples. Furthermore, we conducted differential analysis specifically focusing on the FTO cluster grouping.

To analyze the expression of m6A-related genes in all samples, we first obtained m6A-related genes from the literature (Shi et al., 2019; Chen et al., 2019; Xu et al., 2020; Du et al., 2019), including 11 writer genes, 23 reader genes, and 3 eraser genes, totaling 37 genes. We then intersected the existing expression profiles for 27 genes. First, we used the R package “heatmap” to draw the expression heat map of the genes in all samples, and then we used the “ggpubr” package to draw grouped boxplots based on the normal and patient samples. The Wilcoxon rank sum test method was used for statistical significance between groups, and p < 0.05 was considered statistically significant. The “RCircos” package was used to draw a location map of the 27 genes on the chromosome (Zhang et al., 2013). Chromosomal data were provided by the R package, and the location information of the genes on the chromosome was downloaded from the Ensembl database (Yates et al., 2000).

To further analyze the correlation between writer and eraser gene expression in all patients, Pearson’s correlation was calculated between the two genes. The absolute value of the correlation coefficient was greater than 0.5 and the p value was less than 0.05. We used the R package “ggplot2” to draw a scatterplot of the correlations between gene pairs that met the requirements and fit the correlation curve and used the “ggExtra” package to draw a histogram of the edges of the graph.

Owing to the important influence of the m6A modification process, there may be different m6A modification states between normal and patient samples; therefore, we constructed two training groups based on m6A-related genes from GSE30119 as diagnostic models. In addition, we utilized GSE16129 as an external validation dataset. Here, we screened all m6A genes via least absolute shrinkage and selection operator (LASSO) regression using the R package “glmnet” and selected the best lambda value. Only genes with non-zero coefficients were retained. The genes used to construct the model, and their corresponding coefficients were displayed in the form of forest plots using the R package “forestplot.” Subsequently, to reveal the common effect of m6A gene expression on the diagnostic performance of the model, the R package “rms” was used to construct a logistic regression model of the m6A genes with the most significant weights in the previous LASSO model and visualized using a nomogram. To verify the predictive power of the diagnostic model, the R package “pROC” was used to draw single-gene receiver operating characteristic (ROC) curves and calculate the area under the curve (AUC) (Robin et al., 2011). An internal dataset and decision curve analysis (DCA) were used to illustrate the validity of the nomogram; the curves were drawn using the “ggDCA” R package.

The expression of different genes, especially those that regulate the same biological processes, is interrelated. To reveal the relationship between m6A-related genes, a PPI network was constructed based on m6A-related genes using the STRING database with the above genes as inputs and the confidence threshold at the default value of 0.4 (Von Mering et al., 2003). Subsequently, the PPI network was exported and further analyzed using Cytoscape software (Shannon et al., 2003). The network attributes of each node were calculated, and the plug-in cytoHubba was used to mine the hub nodes based on their degree (Chin et al., 2014). The 10 nodes with a degree of TOP10 were defined as hub nodes. These nodes exhibited a high level of connection with other nodes; therefore, they may play an extremely important role in the regulation of the entire biological process and warrant further study. Therefore, we conducted further prediction research on the 10 hub nodes based on the miRNet database to predict the miRNAs and transcription factors of the hub nodes (Chang et al., 2020). The predicted results were processed and plotted using Cytoscape software.

All hub genes were screened using LASSO regression, and a LASSO model was constructed and displayed as a forest graph. To verify the predictive power of the diagnostic model, the R package “pROC” was used to draw the ROC curve of the model and calculate the AUC.

Because of the pervasive heterogeneity among samples, unsupervised clustering of samples based on m6A regulators was applied to resolve heterogeneity and reclassify samples. Different m6A modification patterns were identified based on the expression of m6A regulators. Clustering was performed using the R package “ConsensusClusterPlus,” and the number of clusters was estimated (Wilkerson and Hayes, 2010). The basic principle of consensus clustering assumes that samples extracted from different subclasses of the original dataset constitute a new dataset, and that different samples from the same subclass are extracted and clustered on the new dataset. Accordingly, both the numbers of clusters and samples within the class should be similar to those of the original dataset. Therefore, the more stable the resulting cluster is with respect to the sampling variation, the more representative the cluster is of a true subclass structure. The resampling method could disrupt the original dataset, so clustering analysis was performed on each resampling sample, and then the results of multiple clustering analyses were comprehensively evaluated to determine the consensus.

To reveal the biological differences between the two sample groups, we conducted Gene Set Enrichment Analysis (GSEA) and displayed the results using volcano and heat maps. Significant differentially expressed genes (DEGs) were defined as having a corrected p value <0.05 and |log2FC| > 0.5.

We performed Gene Ontology (GO) enrichment analysis on significant DEGs to annotate their functions, focusing on three categories: biological process, molecular function, and cellular component (Ashburner et al., 2000). The Kyoto Encyclopedia of Genes and Genomes (KEGG) database provides valuable information on genomes and pathways (Kanehisa and Goto, 2000). The R package GOPlot and KEGGPlot were used to annotate the GO and KEGG functions of all significant DEGs to identify the enriched biological processes (Yu et al., 2012). GSEA helps determine statistical differences in predefined gene sets between two biological states and assesses changes in pathway activity (Subramanian et al., 2005). We downloaded reference gene sets “c5.go.v7.4.Entrez.GMT” and “c2.cp.kegg.v7.4.Entrez.GMT” from the MSigDB database and performed GSEA using “clusterProfiler,” with p < 0.05 as the significance criterion (Liberzon et al., 2015). Gene Set Variation Analysis (GSVA) is a non-parametric method that converts gene expression matrices into gene set expression matrices to evaluate pathway enrichment across samples (Hänzelmann et al., 2023). To examine biological process variation between the two groups, we calculated enrichment scores for each pathway from the “c2.cp.kegg.v7.4.Entrez.GMT” dataset and used the “limma” package to identify significantly differential pathways. GSVA results were visualized using “pheatmap” considering p < 0.05 as statistically significant.

To further explore the similarities and differences in immune cell infiltration levels between the two groups of samples, the “GSVA” package was used following the single-sample GSEA (ssGSEA) method. The marker genes of the 28 immune cells were obtained from the literature and used as the background gene set for ssGSEA of each sample (Charoentong et al., 2017). The infiltration of all immune cells was visualized using box plots. Simultaneously, the R package “corrplot” was used to draw a correlation map between the immune cells for the two groups of samples to reveal the similarities and differences in the degree of correlation of immune cells in different cancer states. In addition, to directly view the correlation between the hub genes and the level of immune cell infiltration, a correlation scatter plot was drawn for the gene–immune cell pairs with significant correlations, and a correlation curve was fitted.

To maximize the accuracy of the results, the R package “CIBERSORT” was used to evaluate the infiltration level of immune cells (Steen et al., 2020), and the content of 22 immune cells in each sample was calculated based on the LM22 background gene set provided by CIBERSORT to reflect the infiltration level. CIBERSORT is based on the principle of linear support vector regression to deconvolute the transcriptome expression matrix and estimate the composition and abundance of immune cells in mixed cells. The results were displayed using heat maps and stacked bar charts drawn using the R package “ggplot2.” For gene–immune cell pairs with significant correlations, we drew a correlation scatterplot and fit a correlation curve. Samples with p < 0.05 were included to obtain the immune cell infiltration matrix.

In order to analyze the correlation between these hub genes, a correlation heat map was drawn using the R package “corrplot.” In addition, in order to further study the correlation of these genes with endoplasmic reticulum stress and mitophagy process, relevant genes were retrieved from the GeneCards database with the keywords “endoplasmic reticulum stress” and “mitophagy,” and then the correlation between hub genes and these genes was calculated and visualized in the form of a bubble chart (Stelzer et al., 2016). For the most significant hub gene pairs, a correlation scatter plot was drawn and a correlation curve was fitted.

Staphylococcus aureus, a pathogen causing asteomyelitis, preserved and provided by the laboratory (Guangdong Provincial Key Laboratory of Bone and Cartilage Regeneration Medicine, Nanfang Hospital, Southern Medical University Guangzhou, Guangdong, China). Prior to conducting the infection experiments, S. aureus was added to 10 mL of fresh tryptic soy broth and incubated overnight at 37°C with agitation at 120 rpm. After centrifugation, the bacteria were washed three times with phosphate-buffered saline (PBS) and resuspended in PBS. The concentration of S. aureus was adjusted to an optical density of 0.5 at 600 nm, which is approximately equivalent to 1 × 10⁸ colony-forming units per milliliter (CFU/mL), to ensure consistent inoculum densities (Lu et al., 2023). The resulting bacterial suspension was appropriately diluted for infection experiments in a RAW 264.7 macrophage cell line (Procell, Wuhan, China). Cells were seeded at a density of 1 × 10⁶ cells/well in a 6-well plate and cultured in Dulbecco’s modified Eagle’s medium (PM150210, Procell) supplemented with 10% fetal bovine serum (FBS; 164,210–50, Procell) and 1% penicillin–streptomycin (PB180120, Procell). To evaluate the gene expression response of macrophages to S. aureus infection, cells were infected with S. aureus at doses of 100, 10, and 1× the multiplicity of infection (MOI). After 1 h of infection, cells were treated with 20 μg/mL gentamicin (215–778-9, Sigma-Aldrich, St. Louis, MO, United States) for 30 min to kill any remaining extracellular S. aureus, thus eliminating extracellular bacteria. After washing thrice with PBS, the cells were incubated in fresh culture medium containing 10% FBS for an additional 24 h. Subsequently, the RNA was collected for mRNA expression analysis.

RAW 264.7 cells were lysed using TRIzol reagent (TaKaRa Bio, Kusatsu, Japan) to isolate RNA. Subsequently, the RNA was reverse transcribed into cDNA using the Uni All-in-One First-Strand cDNA Synthesis SuperMix for qPCR kit (TransGen Biotech, Beijing, China). qPCR was performed using CFX96 (Bio-Rad, Hercules, CA, United States) and RealStar Power SYBR qPCR Mix (GenStar, Beijing, China). This method allows for the detection and quantification of gene expression levels in cDNA samples. To ensure accurate normalization and comparison of gene expression data, the relative expression of the target genes was normalized to that of the reference gene, GAPDH. Data analysis and quantification were performed using the 2^–ΔΔCt method, which calculates the fold change in gene expression relative to a control sample. The PCR primers were 5′-GACACTTGGCTTCCTTACCTG-3′ (FTO forward) and 5′-CTCACCACGTCCCGAAACAA-3′ (FTO reverse); 5′-AGGTCGGTGTGAACGGATTTG-3′(GAPDH forward) and 5′-GGGGTCGTTGATGGCAACA-3′ (GAPDH reverse).

To collect RAW 264.7 cells infected with S. aureus, Cell Lysis Buffer for Western or IP (P0013, Beyotime, Shanghai, China) supplemented with a protease inhibitor cocktail (ST506-2, Beyotime), was used for protein extraction and purification. Equal amounts of protein were separated by SDS-PAGE on a 10% polyacrylamide gel and then transferred onto a PVDF membrane. The membrane was blocked with skim milk at room temperature for 1 h and washed thrice with TBST. Subsequently, the membranes were incubated overnight with primary antibodies against FTO (1:1000; 41,548, Signalway Antibody LLC, College Park, MD, United States), β-actin (1:1000; 52,901, Signalway Antibody LLC), Anti-IL1B (Signalway Antibody LLC, dilution 1:1000), Anti-IL6 (Signalway Antibody LLC, dilution 1:1000), Anti-NFKB (Signalway Antibody LLC, dilution 1:1000), Anti-p-NFKB (Signalway Antibody LLC, dilution 1:1000), and Anti-FOXO1 (Signalway Antibody LLC, dilution 1:1000). The following day, the membrane was washed thrice with TBST and incubated with horseradish peroxidase-conjugated goat anti-rabbit IgG secondary antibody (1:1000; L3012, Signalway Antibody LLC). Finally, images were captured using GelView 6000Plus (Biotend, Guangzhou, China).

All data processing and analysis was performed using Excel (Microsoft, Redmond, WA, United States) and R software (version 4.0.2). For comparisons of two groups of continuous variables, the statistical significance of normally distributed variables was estimated using an independent Student’s t-test, and the differences between non-normally distributed variables were analyzed using the Mann–Whitney U test (i.e., Wilcoxon rank sum test). The chi-squares test or Fisher’s exact test was used to compare and analyze the statistical significance of categorical variables between the two groups. The Kruskal–Wallis test was used to compare two or more groups, and the Wilcoxon test was used to compare two groups. ROC curves were drawn using the pROC package in R, and the AUC was calculated to assess the accuracy of the risk score in estimating prognosis. All statistical p values were two-sided, and p < 0.05 considered statistically significant. Two-tailed p values <0.05 was considered statistically significant.

We first examined the gene expression distribution of the original expression profile of the GSE30119 dataset before and after batch effect correction. The samples showed a serious batch effect when integrated directly, and samples from different sources showed significantly different expression distribution characteristics. After batch effect correction and log normalization, the expression distribution of all samples tended to be consistent, improving the accuracy and robustness of the downstream analysis (Figure 2A). The data spectrum expression distribution of the GSE30119 dataset before and after standardization correction is shown in the box diagram. The heatmap and volcano map show the difference in expression between normal and diseased samples in the dataset GSE30119 (Figures 2B,C).

The expression heterogeneity of all m6A genes in normal and diseased samples was visualized using a heat map and a grouping box diagram (Figures 2E,F). The results showed that among the three m6A gene types, the difference between the readers in the two groups was more obvious than that between writers and erasers. To construct a panorama of m6A-related genes in all samples, we examined the localization of these genes on chromosomes and found that some genes were very close to each other, indicating that they were closely related at the genomic level (Figure 2D). The results showed that some genes were very close to each other on the chromosome, indicating that these genes were closely related at the genomic level. To further analyze the relationship between the expression of m6A writer and eraser genes, we calculated the correlation between these genes and obtained the correlation coefficients R > 0.5 and p < 0.05, indicating positive correlations (Supplementary Figure S1).

Since the expression of m6A-regulated genes has important biological significance, we constructed a diagnostic model for S. aureus infection based on all m6A genes.

First, the 27 m6A genes in the training and validation sets were regressively screened using LASSO, and the best lambda values were obtained for 17 genes in training set C1 (Figures 3A,B) and the 7 genes in validation set C2 (Figures 3D,E). Subsequently, forest maps were used to visualize the effects of these m6A gene-containing diagnostic models in the training (Figure 3C) and validation sets (Figure 3F).

Subsequently, to verify the accuracy of the model, we created nomograms and found that in the GSE30119 training set C1, the patient prediction risk score correlated with the disease risk of patients with S. aureus infection, highlighting the accuracy of the model (Supplementary Figure S2A). The predictive performance of the model was further validated using recall curves and decision curve analysis (DCA), which confirmed that the model exhibited superior predictive performance and robustness (Supplementary Figures S2C,E). Additionally, the nomogram model, recall curves, and DCA demonstrated the same efficacy as observed in the C1 validation set (Supplementary Figures S2B,D,F).

To further identify the key m6A regulators for diagnosing patients with the disease, we compared the factors with diagnostic significance in the two sets and found that only FTO was differentially expressed (Figures 4A,B). The AUC for FTO as a diagnostic marker for S. aureus infection was 0.857 in the training set and 0.886 in the validation set (Figures 4E,F). Moreover, the difference in FTO expression between different disease groups was significant, suggesting that FTO has high diagnostic value. Combined with the above forest diagram, it was confirmed that FTO is a low-risk gene for S. aureus infection.

To validate and analyze whether the diagnostic models based on the LASSO algorithm for C1 and C2 could independently distinguish between diseases, we plotted ROC curves and calculated the AUC. The AUC for C1 was 0.880, and for C2, it was 0.861, indicating good diagnostic efficacy (Figures 4C,D). Furthermore, to enhance the validation of diagnostic performance, we selected the GSE16129 dataset related to S. aureus infection osteomyelitis for external validation. LASSO regression, forest plots, and nomograms all confirmed that FTO demonstrated diagnostic efficacy (Supplementary Figures S3A–D). Additionally, calibration curves and decision curve analysis (DCA) further supported these findings (Supplementary Figures S3E,F).

To investigate the role of the key diagnostic marker FTO in disease, we conducted unsupervised consensus clustering on 143 samples from the GSE30119 dataset based on gene expression. The disease samples were clustered according to different m6A regulator subclasses, leading to an optimal k value of 2 (Figures 5A–C), which allowed us to classify the samples into two distinct subtypes (A: n = 70; B: n = 73, Figure 5D). We then compared the cluster groupings with high- and low-expression FTO levels, visualizing the results in a Sankey diagram (Figure 5F). Notably, Class A predominantly represented the disease samples, suggesting that FTO may play a significant role in differentiating between FTO Cluster A and FTO Cluster B.

To further analyze the influence of FTO on diseases, we analyzed the differences between groups with high and low FTO expression and obtained a total of 645 co-expressed genes that may be subject to the same regulatory processes and reflect similar biological functions (p value <0.05 and | log2FC | > 1) (Figure 5G). Additionally, I selected the top 50 genes for visualization in a heatmap (Figure 5E).

To gain a deeper understanding of FTO-related differential genes and their biological significance, we utilized the STRING database was used to construct a protein–protein interaction network of differential genes for the 645 FTO co-expressed genes (Figure 5H). The tightly linked genes of the PPI network module were screened using the MCODE plug-in in CYTOSCAPE (v3.7.2), and the highest confidence interaction score was set to 0.4 and visualized (Figure 5I). Simultaneously, the cytoHubba plug-in was used to screen the Top 30 closely linked genes (Figure 5J). The intersection of the two methods is shown in a Venn diagram, which shows 19 closely related FTO co-expressed genes (Figure 5K). These genes are closely related to FTO expression.

Subsequently, to further understand the interaction of the 19 FTO co-expressed genes at the post-transcriptional stage, we identified 589 miRNAs and 206 transcription factor regulatory genes targeted by them and constructed a network (Figure 5L).

GO analysis showed enrichment in DNA binding, mismatch repair complex binding, and cell cycle. Function-related genes were enriched in immune cell proliferation, differentiation, activation, and oxidation (Supplementary Figure S4). The KEGG enrichment results indicated that the entries were enriched in various cellular senescence, tumor, and immune cell-related pathways. GSEA enrichment analysis was performed between the groups with high and low FTO expression (c2.all.v7.2.symbols.gmt was used as the background set) (Supplementary Tables S1–S3). The results suggested that the related molecular mechanisms mediated by PD-1, TH1TH2, CTLA4, and other pathways were significantly enriched in patients with high FTO expression (Figure 6A).

To further investigate the biological characteristics of FTO expression in diseased tissues, we performed unsupervised consensus clustering using the expression of 19 hub co-expression molecules associated with FTO. Cluster analysis was conducted on disease sample data to classify the samples into different subclasses based on these molecules. After evaluating multiple cluster analyses for consistency, we selected an optimal k value of 3 based on the delta plot, resulting in three distinct subtypes (A: n = 26; B: n = 58; C: n = 14; Figures 6B–D).

Subsequently, we performed Gene Set Variation Analysis (GSVA) between the groups, focusing on comparisons among Cluster C vs. Cluster A, Cluster C vs. Cluster B, and Cluster B vs. Cluster A, using the c2.all.v7.2.symbols.gmt background set. The differentially enriched pathways identified included amino acid metabolism and immune- and infection-related processes (Figures 6E–G).

Immune cell infiltration was assessed using the GSE30119 dataset. The CIBERSORT algorithm is shown in a box chart to evaluate the difference in immune cell infiltration between the normal and disease sample groups (Figure 7A), and the immune infiltration correlation analysis was further performed on the disease-related diagnostic prediction model according to the difference in immune cell expression between the FTO high and low expression groups (Figure 7B). The expression of naïve B cells, CD8+ T cells, naïve CD4+ T cells, activated memory CD4+ T cells, gamma delta T cells, resting natural killer cells, monocytes, M2 macrophages, resting mast cells, and neutrophils differed between the normal and disease sample groups. The number of CD8+ T cells, naïve CD4+ T cells, T cell helper cells, and neutrophils differed between the high and low FTO expression groups.

We then performed ssGSEA to estimate the number of specific infiltrating immune cells and the specific immune response activity, which defined an enrichment score to determine the absolute enrichment of a gene set in each sample in the dataset. We compared the differences in immune cell infiltration among the three clusters from FTO hub co-expression molecular typing (Figure 7C) and between the high and low FTO expression groups (Figure 7D).

To analyze the correlation between FTO expression and immune cells in diseased samples, we drew a scatter plot showing the expression correlation scatter plot between the genes with the strongest correlation (correlation coefficient R > 0.6 and significant sex p < 0.01). Spearman correlation analysis was used to determine the correlation between FTO expression, the proportion of immune cells, and immune reactivity. The results showed that FTO expression was closely related to Activated.B.cell (R = 0.548), Immature.B.cell (R = 0.479), T.follicular.helper.cell (R = 0.444), Activated.CD8.T.cell (R = 0.437), MDSC (R = 0.418), and type.17.T.helper cells (R = 0.402), with a positive correlation (Supplementary Figures S5A–F).

To investigate the effect of S. aureus infection on FTO expression in macrophages, RAW 264.7 cells were infected with different MOIs of S. aureus. We observed a significant decrease in both the mRNA and protein levels of FTO as the MOI increased. Following infection with S. aureus at an MOI of 10, the FTO expression in macrophages was statistically significant (Figure 8A). Furthermore, Western blot experiments confirmed that FTO protein expression was reduced in macrophages after infection with S. aureus at an MOI of 10, which was consistent with the observed changes in mRNA expression (Figure 8C). Based on the groups—normal group, S. aureus-infected macrophage group, si_Fto group, and S. aureus-infected macrophage + si_FTO group, we conducted Western blot analysis to assess the expression of inflammatory proteins and the FOXO1/NFKB signaling pathway proteins. We found that in the S. aureus-infected macrophage group and the S. aureus-infected macrophage + si_FTO group, IL6 and IL1β inflammatory protein levels increased following si_FTO treatment. Additionally, we detected the expression of p-NFKB, NFKB, and FOXO1 proteins, suggesting that FTO may regulate the macrophage phenotype through the FOXO1/NFKB pathway (Figure 8D). Additionally, we performed statistical analysis on the protein bands using bar charts, which demonstrated significant differences among the groups (Figure 8B).