Identification of DEGs
In our study, 535 DEGs were identified between SONFH samples and control samples. Among them, 299 were upregulated and 236 were downregulated (SONFH vs. control). The volcano plot and heatmap of gene expression are shown in Figures 2A and 2B.
Construction of co-expression networks
The sample clustering tree indicated that there was no abnormal sample (Figure 3A-B). After calculation, the best soft-thresholding power was set at 18 (Figure 3C). Finally, each module was assigned a color, and a total of 9 modules in GSE123568 (Figure 3D) were identified in this study. Furthermore, the result of the module-feature relationship revealed that the brown module had the highest correlations with SONFH (cor = 0.68, p = 1e−05, Figure 3E). Thus, 850 genes in the brown module were selected for further exploration.
Identification of IRGs-SONFH and functional enrichment analysis
Then, we took the intersection of DEGs, genes in key modules, and IRGs and identified 25 IRGs-SONFH (Figure 4A). To explore the function of 25 IRGs-SONFH in SONFH, the GO terms are shown in Figure 3. In BP analysis (Figure 4B), IRGs-SONFH mainly participated in response to molecule of bacterial origin, neutrophil activation, response to lipopolysaccharide, cellular response to biotic stimulus, and cellular response to molecule of bacterial origin. In CC analysis (Figure 4C), IRGs-SONFH significantly participated in the membrane microdomain, membrane raft, secretory granule membrane, endocytic vesicle, and phagocytic vesicle. MF analysis showed that IRGs-SONFH significantly participated in amide binding, peptide binding, immune receptor activity, pattern recognition receptor activity, and lipopolysaccharide binding (Figure 4D). KEGG analysis was performed to explore the pathways of these 25 IRGs-SONFH. The KEGG terms of IRGs-SONFH are shown in Figure 4E. As shown, these IRGs-SONFH were mainly enriched in lipid and atherosclerosis, tuberculosis, neutrophil extracellular trap formation, TLR signaling pathway, and legionellosis.
Identification of hub genes
The PPI network between IRGs-SONFH was established using the STRING database, interactions of 25 IRGs-SONFH were displayed in Figure 5A. 4 hub genes (CD14, CYBB, NOD2, and TLR1) were identified by MCODE plug-in Cytoscape (Figure 5B).
The correlation analysis between hub genes and the functional similarity analysis of hub genes
The correlation between these 4 hub genes was investigated using the corrplot package, CD14 and TLR1 had the strongest correlation (r = 0.85) (Figure 5C). We analyzed the functional similarity of these hub genes by the “GOSemSim” package in R. The results showed that 3 hub genes, including CD14, NOD2, and TLR1 (similarity score>0.5), had higher functional similarity (Figure 5D).
Validation and efficacy evaluation of hub genes
We explored the expressions of these genes between SONFH and control samples in GSE123568 and found that those genes exhibited higher expression levels in SONFH (Figure 6A). In addition, the relative expressions of the above four hub genes were investigated by qRT-PCR. As shown in Figure 6C-F, the relative expressions of CD14, CYBB, NOD2, and TLR1 were also significantly increased in peripheral blood of SONFH patients compared to controls. Furthermore, we executed a ROC curve analysis to calculate their sensitivity and specificity for the diagnosis of SONFH (Figure 6B). The AUC values of CD14, CYBB, NOD2, and TLR1 were 0.847, 0.753, 0.767, 0.847, respectively, demonstrating that these genes have high sensitivity and specificity for SONFH diagnosis.
Correlation analysis of hub genes and immune cells
To further understand the role of these genes in immune infiltration, we used spearman correlation analysis to determine whether these hub genes were related to immune cell infiltration. Correlation analysis showed that 4 hub genes including CD14, CYBB, NOD2, and TLR1 had significantly positive relationship with Type 1 T helper cell, T follicular helper cell, regulatory T cell, plasmacytoid dendritic cell, neutrophil, natural killer T cell, natural killer cell, monocyte, memory B cell, myeloid-derived suppressor cell, mast cell, macrophage, immature dendritic cell, immature B cell, gamma delta T cell, eosinophil, effector memory CD8 T cell, effector memory CD4 T cell, central memory CD4 T cell, central memory CD8 T cell, activated dendritic cell (Figure 7).
The function of our hub genes was explored via GSEA. Genes in the high-expression cohorts of CD14, and TLR1were all highly enriched in leishmania infection, Toll-like receptor signaling pathways and Fc gamma R-mediated phagocytosis (Figure 8A, 8D). Genes in the high-expression cohorts of CYBB, and NOD2 were all highly enriched in spliceosome, lysosome and B-cell receptor signaling pathways (Figure 8B, 8C). Genes in the low-expression cohorts of CD14, CYBB, NOD2, and TLR1 were all highly enriched in olfactory transduction, linoleic acid metabolism, and basal cell carcinoma (Figure 8). After considering the results of GSEA, we concluded that these four genes might be highly correlated with immune and inflammation.
A total of 17 potential drugs for treating SONFH patients were identified when the drug-gene interactions were explored using DGIdb (Table 2). Additionally, drug-gene networks were constructed by Cytoscape (Figure 9A). However, we did not find any small molecule drugs that could target TLR1 in this database.
Prediction of key miRNAs and TF
The miRNA and TFs regulatory network of the 4 hub genes was constructed using miRNet. As illustrated in Figure 9B, the interaction network consisted of 4 hub genes and 59 miRNAs. Specifically, 9 miRNAs (ie, hsa-mir-335-5p, hsa-mir-100-5p, hsa-mir-3687) targeting CD14, 28 miRNAs (ie, hsa-mir-6826-3p, hsa-mir-6845-3p, hsa-mir-6859-3p) targeting CYBB, 12 miRNAs (ie, hsa-mir-215-5p, hsa-mir-122-5p, hsa-mir-320a) targeting NOD2, 15 miRNAs (ie, hsa-mir-34a-5p, hsa-mir-3662, hsa-mir-4511) targeting TLR1. The interaction network consisted of 4 hub genes and 30 TFs. We found that 12 TFs (ie, CEBPB, FOS, JUN) could regulate CD14. 8 TFs (ie, NFIC, NFYA, YY1) could regulate CYBB. 11 TFs (ie, MAX, USF1, USF2) could regulate NOD2. 8 TFs (ie, MEF2A, HINFP, TP63) could regulate TLR1.
CeRNA regulatory network construction
To elucidate the potential molecular mechanism of lncRNAs in SONFH, we constructed a lncRNA-miRNA-mRNA interaction network. Briefly, 20 DEmiRNAs with the threshold criterion of adjusted p<0.05 were screened by GSE89587. 1 miRNA (hsa-miR-320a) was obtained by taking the intersection of 64 miRNAs predicted by the hub genes above and 20 DEmiRNAs. We used the database miRNet to predict the lncRNAs that interacted with the selected miRNAs (hsa-miR-320a). Finally, we obtained a ceRNA network which included 67 lncRNAs, 1 miRNA (hsa-miR-320a), and 1 mRNA (NOD2) (Figure 9C).