We obtained the AS datasets (GSE25101, GSE73754) and the LBMD datasets (GSE56815, GSE2208) from the GEO database (Table 1). All analyses were performed using R software (version 4.0.2).

Overlapping genes between each disease dataset were identified. For AS, the intersection of the genes analyzed in GSE25101 and GSE73754 was determined, and the same for LBMD (GSE56815, GSE2208). The expression of these overlapping genes was extracted from each dataset. We used the ComBat method to eliminate batch bias using the sva package (version 3.36.0).

The corrected AS and LBMD datasets were subjected to differential expression analysis (DEA) by the limma package (version 3.44.3). We defined the genes with p-values < 0.05 as differential expressed genes (DEGs). The ggplot2 package (version 3.3.2) and pheatmap package (version 1.0.12) were used to depict gene expression by heatmaps and volcano plots.

The immune-related genes (IRGs) were obtained from the ImmPort database (Bhattacharya et al., 2018). The AS DEGs, LBMD DEGs, and IRGs were taken intersection to obtain the dysregulated IRGs in both AS and LBMD, and we called them immune-DEGs (I-DEGs).

The correlation between I-DEGs and DEGs in AS and LBMD was analyzed by calculating the Pearson coefficient in the Hmisc package (version 4.4.1). The DEGs related to I-DEGs were obtained with moderate correlation (p-values < 0.05 and | r | > 0.40). In addition, we constructed a co-expression network of I-DEG-related DEGs with strong correlation (p-values < 0.05 and | r | > 0.70).

Louvain is a clustering algorithm frequently employed in network clustering, notably in human PPI networks, intending to decrease external connections while encouraging intra-community connections (Smith et al., 2020). We clustered the co-expression network using the Louvain algorithm in the igraph package (version 1.2.5), where the weight was set to | r |. The modules with less than 30 genes were deleted, and the Cytoscape software (version 3.9.0) was used for visualization.

Module eigengenes (MEs) were defined as the genes in the first principal component of gene modules using the PCA algorithm. The I-DEGs were defined as the immune module as a whole. The correlation between the MEs of I-DEG-related DEGs modules and the ME of the immune module was analyzed by the Pearson correlation coefficient and visualized by the corrplot package (version 0.84). Among them, the modules with the highest correlation were considered the key modules and studied in subsequent analyses.

The Metascape database (Zhou et al., 2019) was used to perform GO and KEGG functional enrichment analysis for the key module genes. The criteria of the analyses were p-values < 0.01, min (overlap) = 3, and min (enrichment) = 1.5, and the outcomes were visualized using heatmaps.

The PPI network of the key modules was constructed by the STRING database (Szklarczyk et al., 2021) and visualized by Cytoscape. The cytohubba plug-in was used to find the top10 hub genes of the key modules.

The I-DEGs with the same expression trend in AS and LBMD were considered the key I-DEGs. They were imported into the STRING database with the hub genes and visualized by Cytoscape to explore the relationship between the key I-DEGs and hub genes.

To explore the function of key I-DEGs as a whole in more depth, we designated them as the key signature. The AS and LBMD datasets were split into training and testing sets with a 7:3 ratio. The Glmnet package (version 4.0.2), randomForest package (version 4.6.14), XGBoost package (version 1.4.1.1) and e1071 package (version 1.7.3) were used to construct the diagnostic model with the key signature. We get the LASSO model and the corresponding feature coefficients by performing the least absolute shrinkage and selection operator (LASSO) binomial analysis in the Glmnet package. The penalty parameter (λ) was decided by the minimum criteria. The randomForest package was used to build a random forest model for the key signature, and the decreasing accuracy method (Gini coefficient method) was used to obtain the feature importance values. The XGBoost package was used to build the XGBoost model. The XGBoost model was a tree-based algorithm that can provide each feature importance score and rank them. The e1071 package was used to build the support vector machine (SVM) model, a typical supervised machine learning method.

To explore the interaction context of key I-DEGs, especially in the skeletal system, we imported key I-DEGs into the IID database and constructed the PPI network in all tissue and bone. We called the genes that connected key I-DEGs bridge genes, which may help us understand the context of key I-DEGs and make the gene interaction connection abundant.

To investigate the pathogenic mechanism and better comprehend the function of key I-DEGs, the key I-DEGs were imported into the GO (Carbon et al., 2009) and KEGG database (Kanehisa and Goto, 2000) to obtain human biological processes and pathways comprising at least two genes. To provide a better explanation for the findings, the biological processes acquired were combined and clustered by ClueGo (Bindea et al., 2009). KEGG database was utilized to explore the common pathways map, getting potential key I-DEG and pathway associations.

The CIBERSORT (Newman et al., 2015) algorithm was used to acquire the immune infiltration matrix in the AS datasets. Using the ggplot2 package, each sample and group’s immune infiltration was shown graphically. The differences between the two groups were examined by the Wilcoxon test. Using the corrplot package, we plotted correlation heatmaps to demonstrate the association between the 22 immune cells and between key I-DEGs and immune cells.

After batch correction, the differences between samples in AS and LBMD datasets were significantly reduced (Figure 2). Based on DEA, we got DEGs of AS and LBMD, including 1,492 up-regulated and 1818 down-regulated genes in AS, and 730 up-regulated and 454 down-regulated genes in LBMD. Figure 3 showed the expression of DEGs. From the ImmPort database, we obtained 1793 IRGs (Supplementary Table S1). Figure 4 showed that a total of 19 I-DEGs were identified. They were not only common DEGs of AS and LBMD, but also IRGs. These I-DEGs were used for our subsequent analysis in depth.

Through the Pearson correlation analysis, we obtained 2,227 DEGs related to I-DEGs in AS and 735 in LBMD (Supplementary Table S2). Furthermore, the co-expression networks of these genes were constructed in AS and LBMD, respectively. By the Louvain community discovery algorithm, the gene modules with less than 30 genes were deleted, and the gene co-expression network of AS was re-clustered into seven modules. Similarly, LBMD was re-clustered into three modules (Figures 5A,B). Density is an assessment metric for measuring the connectivity degrees of the network. The density of these modules was shown in Table 2. It was noted that the density of these modules was higher than that of the co-expression network (AS: 0.01673044, LBMD: 0.05024927). It demonstrated that the multilevel algorithm’s partitioning result was reliable. Figures 5C,D showed that M4 was the most related to the immune module in AS and M3 in LBMD. These modules were called key modules and were closely related to immunity which may play an essential role in AS and LBMD.

Figure 6A showed that the genes of AS key module AS-M4 were mainly involved in several GO BPs, for example, cellular carbohydrate metabolic process, leukocyte activation, secretion, protein phosphorylation, positive regulation of cell death, myeloid leukocyte activation, and inflammatory response. Meanwhile, AS-M4 genes were mainly enriched in several pathways (Figure 6B), Neutrophil extracellular trap formation, Lipid and atherosclerosis, Necroptosis, Yersinia infection, Pathways in cancer, Non-small cell lung cancer, Legionellosis, Osteoclast differentiation, Proteoglycans in cancer, Apoptosis. This module was closely related to immunity, such as immune cell activation and inflammatory response. Interestingly, osteoclast differentiation was one of its significant pathways, confirming its vital role in affecting bone homeostasis.

Figure 6C showed that the biological processes of LBMD key module LBMD-M3 were mainly enriched in mRNA processing, ncRNA processing, regulation of cellular amide metabolic process, regulation of mRNA metabolic process, organelle fusion, cellular response to organonitrogen compound, translation, regulation of I-kappaB kinase/NF-kappaB signaling, RNA 3′-end processing, regulation of circadian rhythm, cell division, protein localization to organelle, cellular response to chemical stress, etc. Meanwhile, LBMD-M3 genes were mainly enriched in several pathways (Figure 6D), mRNA surveillance pathway, Insulin resistance, Ubiquitin mediated proteolysis and Spliceosome. Interestingly, we found this module participated in the regulation of I-kappaB kinase/NF-kappaB signaling, while the receptor activator of NF-kappaB (RANK) had been fully recognized as causing the osteoclast precursor to differentiate into a preosteoclast, which was an essential mechanism of bone destruction (Weitzmann, 2017).

The PPI network of AS-M4 was obtained from the STRING database, including 289 nodes and 787 edges, and the PPI network of LBMD-M3, including 56 nodes and 89 edges (Supplementary Figures S1A,C). The top10 hub genes were identified by the cytohubba plug-in (Supplementary Figures S1B,D).

We further explored the expression of I-DEGs and found that five genes (TNF, CCL3, IFNAR1, PIK3CG, and PTGER2) showed the same trend in AS and LBMD. By examining the expression of IRGs in AS and LBMD, we found that CCL3, PTGER2, and TNF showed a common down-regulation trend in both diseases, while IFNAR1 and PIK3CG showed a common up-regulation trend. These five genes were considered key I-DEGs.

By building the PPI networks of key I-DEGs and hub genes, respectively (Figure 7), we found that key I-DEGs had a good relationship with AS-M4 hub genes, while TNF can be related to UBE2N and TIA1 in LBMD.

In AS, the LASSO model composed of five features was constructed according to the optimal λ value (LASSO score: 1.92exp (CCL3) + 2.42exp (IFNAR1) + 1.25exp (PIK3CG) + 0.2exp (PTGER2)-4.84exp (TNF) (Figures 8A,B). In LBMD, four features constituted the LASSO model (LASSO score: 0.52exp (IFNAR1) + 1.93exp (PIK3CG)-0.65exp (PTGER2)-2.31exp (TNF)) (Figures 8F,G).

We found that key I-DEGs as a whole had a good diagnostic value in LBMD, among which the SVM model had the highest AUC value of 0.76 (Figure 8H). Interestingly, the signature also had an excellent diagnostic value in AS (Figure 8C), further illustrating the critical relationship between AS and LBMD. The specific AUC values of the four models were shown in Table 3.

In LBMD RF and XGBoost models, we found that both TNF and PIK3CG were the top2 most important feature (Figures 8D,E), and PIK3CG was the most important in AS RF and XGBOOST models, and TNF was the second in RF and the fourth in XGBoost (Figures 8I,J). It undoubtedly confirmed the role of TNF and PIK3CG in AS and LBMD, and was worthy of in-depth discussion.

We identified that key I-DEGs potentially can interact with one another through bridge genes (Figure 9A). Furthermore, it was found that IFNAR1, PIK3CG, TNF, and PTGER2 can interact in the skeletal system (Figure 9B). This complements the interaction between key I-DEGs, especially in bone homeostasis.

To get a deeper understanding of the function of key I-DEGs, we integrated them into the GO database and identified 305 biological processes, including at least two key I-DEGs (Supplementary Table S3). The results of ClueGO indicated that 305 biological processes were fused into 131 based on the relationship between GO terms, with MAPK cascade and positive control of acute inflammatory response comprising the majority (Figure 10A).

To further investigate the specific mechanism of key I-DEGs, the key I-DEGs were imported into the KEGG database, and we identified 24 pathways comprising at least two genes after filtering. (Figure 10B). Interestingly, TNF and IFNAR1 were involved in osteoclast differentiation, suggesting that key I-DEGs may cause bone mass loss in AS.

As the most significant pathway of the AS key module, neutrophil extracellular traps (NETs) formation attracted our attention. By viewing the map of NETs formation pathway in the KEGG database, we found that there were four pathways, including PI3K-Akt signaling pathway, NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, and MAPK signaling pathway, participating in NETs formation, and also closely related to the four key I-DEGs (PIK3CG, IFNAR1, CCL3, and TNF) (Figure 10C). NETs formation was not only involved in the pathogenesis of AS, but may also be the key to these key I-DEGs connections.

We analyzed the immune infiltration of AS datasets using the CIBERSORT algorithm to further investigate immune cells’ function in AS and the impact of key I-DEGs on it. Figure 11A provided immune infiltration of the samples. Five immune cell subpopulations revealed significant differences between case and control groups (p-values < 0.05). Compared with the control group, the case group contained a higher proportion of neutrophils, lower T cells CD8, T cells gamma delta, NK cells activated and mast cells activated, and neutrophils were the most significant (p-value < 0.001) (Figure 11B). Figure 11C showed that all the key I-DEGs negatively correlated with macrophage M2, and IFNAR1, PIK3CG, PTGER2 and TNF were positively correlated with macrophage M1. It was a rather interesting outcome indicating key I-DEGs can cause an imbalance of macrophages M1/M2 cells in AS, which has been shown to be associated with bone loss (Wang et al., 2022). It suggests that key I-DEGs can affect bone remodeling by regulating macrophage polarization in AS. From the correlation between immune cells, we found a negative correlation between Macrophages M1 and Macrophages M2, which confirms the imbalance of M1/M2 on the other hand (Supplementary Figure S2). Meanwhile, CCL3, PIK3CG, PTGER2 and TNF negatively correlated with neutrophils. Considering that NETs formation was the most significant pathway in AS key module, neutrophils seem to become more remarkable.

This is the first research to use bioinformatics to evaluate the function of immune-related genes in AS-induced LBMD and investigate the mechanism. Innovative methods such as community discovery and machine learning make the research more exhaustive. We created a comprehensive combination of key I-DEGs and their associated mechanism, filling the gap of the mechanism in previous research and generating novel concepts for future research. Nonetheless, this work was essentially performed by computer technology such as bioinformatics and machine learning, and further experimental validation of this hypothesis is required.