Screening of differentially expressed genes
A total of2036 DE-SONFHGs between the SONFH samples and control samples were identified, and these DE-SONFHGs included 1383 upregulated and 653 downregulated genes. The distribution of these DE-SONFHGs is presented Figure 1.
WGCNA
The overall clustering of the dataset samples was good, so no samples were removed (Figure 2a). Then, the traits of the samples were classified, and sample clusters and clinical feature heatmaps were generated (Figure 2b). According to the position of the blue line, the power threshold was determined to be 24 (R^2 = 0.85), and 48 modules were obtained (Figure 2c-f). Subsequently, similar modules were analyzed and merged with the dynamic cutting tree algorithm with MEDissThres equal to 0.2, and 16 modules remained after merging (Figure 2g). The module with the strongest positive correlation (bisque4 module) and the module with the strongest negative correlation (darkgreen module) were selected as key modules; the bisque4 module included 2382 genes, and the darkgreen module included 245 genes (Figure 2g). Furthermore, 896 hub genes were authenticated in the bisque4 module, and 145 hub genes were authenticated in the darkgreen module (Figure 2h). Finally, a grand total of 1041 hub genes were identified. Moreover, 106 DE-CMR-SONFHGs were identified by taking the intersection of 1041 hub genes, 2036 DE-SONFHGs and 2062 CMRGs (Figure 2i).
Functional enrichment analysis of DE-CMR-ONFHGs
The 106 DE-CMR-SONFHGs were enriched in 143 GO-Biological Process (GO-BP), 50 GO-Cellular Component (GO-CC), 34 GO-Molecular Function (GO-MF) terms and 51 KEGG pathways (Figure 3a and b). These genes were mainly enriched in various GO terms related to protein binding, membrane and stimulus response, such as identical protein binding, plasma membrane and inflammatory response. Furthermore, these DE-CMR-SONFHGs were mainly enriched in phagosomes, the PI3K-Akt signaling pathway, and neutrophil extracellular trap formation.
The linkage between proteins
A PPI network was constructed to explore the interactions of 106 DE-CMR-SONFHGs using the STRING website, 387 protein interaction pairs, including 92 nodes, were identified (Figure 4). TLR4 had the greatest connectivity in the PPI network. In addition, 15 downregulated genes were included in the PPI network.
Screening for key genes by machine learning algorithm
LASSO logistic regression and the SVM algorithm were utilized to identify key DE-CMR-SONFHGs. When lambda.min was 0.1073883, the error rate was the lowest, and 2 candidate genes, namely, PNP and SLC2A1, were identified by LASSO (Figure 5a and b). Additionally, 8 candidate genes, including RNASET2, PNP, SLC2A1, REXO2, CYBA, SOAT1, TFDP1 and LYZ, were selected by SVM (Figure 5c and Table 1). Finally, 2 genes (PNP and SLC2A1) were identified in the results of both LASSO and SVM (Figure 5d). The ROC curves revealed that the area under the curve (AUC) values of PNP and SLC2A1 in the GSE123568 dataset were both above 0.9, indicating that each key gene had diagnostic value in distinguishing between SONFH samples and control samples (Figure 6).
Table 1 Eight candidate genes were screened by SVM, including RNASET2, PNP, SLC2A1, REXO2, CYBA, SOAT1, TFDP1 and LYZ.
FeatureName
|
FeatureID
|
AvgRank
|
RNASET2
|
80
|
6.6
|
PNP
|
67
|
7
|
SLC2A1
|
87
|
9.6
|
REXO2
|
79
|
10.2
|
CYBA
|
32
|
12
|
SOAT1
|
89
|
12
|
TFDP1
|
19
|
13.4
|
LYZ
|
51
|
14
|
Single gene enrichment analysis
To investigate the underlying functions of PNP and SLC2A1, GSEA was performed. PNP was mainly enriched in ‘positive regulation of smoothened signaling pathway’, ‘positive regulation of cytokine production involved in immune response’, ‘Glycosphingolipid biosynthesis lacto and neolacto series’, and ‘B-cell receptor signaling pathway’ (Figure 7a and b and Tables S1-4). SLC2A1 was mainly enriched in ‘cgmp metabolic process’, ‘vacuolar acidification’, and ‘glycosphingolipid biosynthesis lacto and neolacto series’. (Figure 7c and d and Tables S5-8).
Immune cell infiltration analysis
To explore the abundance of 28 infiltrating immune cell populations in all the samples, ssGSEA was performed between the SONFH and control samples. Twenty immune cell populations, such as MDSCs, eosinophils, and CD56bright NK cells, were significantly different between the SONFH and control samples (Figure 8a). Among these differentially abundant cell populations, the abundance of immature dendritic cells was strongly negatively correlated with the abundance of activated B cells, and the abundance of CD56bright natural killer cells was strongly positively correlated with the abundance of activated B cells (Figure 8b). PNP expression had a significant and the strongest negative association with the abundance of plasmacytoid dendritic cells, and PNP expression had a significant and the strongest positive association with the abundance of CD56bright natural killer cells (Figure 8c). SLC2A1 expression had a significant and the strongest negative association with the abundance of T follicular helper cells, and SLC2A1 expression had the strongest positive association with the abundance of T helper 17 cells (Figure 8d).
Drug prediction analysis
After the drug prediction analysis, 16 drugs in the DGIdb database were predicted to target PNP and SLC2A1; of these drugs, 11 drugs were predicted to target SLC2A1, and 5 drugs were predicted to target PNP (Figure 9).
The relevance of copper death-related genes and differentially expressed genes to disease
To investigate the correlation between copper death-related genes and differentially expressed genes in femoral head necrosis, we calculated correlations based on the expression of the 2 key genes that were identified and the expression of 10 copper death-related genes, and we determined the corresponding p value values and correlation coefficients r. Significantly correlated relationship pairs were identified based on the correlation threshold p value < 0.05 and |r|> 0.3 (Table S9). The correlation results were visualized through the R package ‘ggplot2’ (Figure 10) according to the correlation ranking, where a significant and the strongest positive correlation was observed among SLC2A1, PNP and CDKN2A, and a significant and the strongest negative correlation was observed among SLC2A1, PNP and PDHB.
The expression of the 10 copper death-related genes was extracted from the GSE123568 dataset and combined with the sample grouping information. The expression of copper death-related genes in the disease samples and control samples was evaluated by the wilcox.test method (Figure 11). The results showed that four copper death-related genes (LIPT1, DLD, PDHB, and MTF1) were significantly different and were differentially upregulated in the disease samples.
Verification of the expression of key targets through qRT‑PCR
qRT‒PCR was performed on human blood samples to analyze the expression of PNP and SLC2A1. Compared to the normal group, the expression of PNP (p = 0.0362) and SLC2A1 (p = 0.0094) was notably reduced in the SONFH group (Figure 12).
Compared to the normal group, the expression of LIPT1 (p = 0.4129), DLD (p<0.0001), PDHB (p =0.0149), and MTF1 (p = 0.0326) was notably decreased in the SONFH group (Figures 11-13).