Identification of mitophagy-related biomarkers in human osteoporosis based on a machine learning model

Abstract

Background: Osteoporosis (OP) is a chronic bone metabolic disease and a serious global public health problem. Several studies have shown that mitophagy plays an important role in bone metabolism disorders; however, its role in osteoporosis remains unclear. Methods: The Gene Expression Omnibus (GEO) database was used to download GSE56815, a dataset containing low and high BMD, and differentially expressed genes (DEGs) were analyzed. Mitochondrial autophagy-related genes (MRG) were downloaded from the existing literature, and highly correlated MRG were screened by bioinformatics methods. The results from both were taken as differentially expressed (DE)-MRG, and Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed. Protein-protein interaction network (PPI) analysis, support vector machine recursive feature elimination (SVM-RFE), and Boruta method were used to identify DE-MRG. A receiver operating characteristic curve (ROC) was drawn, a nomogram model was constructed to determine its diagnostic value, and a variety of bioinformatics methods were used to verify the relationship between these related genes and OP, including GO and KEGG analysis, IP pathway analysis, and single-sample Gene Set Enrichment Analysis (ssGSEA). In addition, a hub gene-related network was constructed and potential drugs for the treatment of OP were predicted. Finally, the specific genes were verified by real-time quantitative polymerase chain reaction (RT-qPCR). Results: In total, 548 DEGs were identified in the GSE56815 dataset. The weighted gene co-expression network analysis(WGCNA) identified 2291 key module genes, and 91 DE-MRG were obtained by combining the two. The PPI network revealed that the target gene for AKT1 interacted with most proteins. Three MRG (NELFB, SFSWAP, and MAP3K3) were identified as hub genes, with areas under the curve (AUC) 0.75, 0.71, and 0.70, respectively. The nomogram model has high diagnostic value. GO and KEGG analysis showed that ribosome pathway and cellular ribosome pathway may be the pathways regulating the progression of OP. IPA showed that MAP3K3 was associated with six pathways, including GNRH Signaling. The ssGSEA indicated that NELFB was highly correlated with iDCs (cor = -0.390, p < 0.001). The regulatory network showed a complex relationship between miRNA, transcription factor(TF) and hub genes. In addition, 4 drugs such as vinclozolin were predicted to be potential therapeutic drugs for OP. In RT-qPCR verification, the hub gene NELFB was consistent with the results of bioinformatics analysis. Conclusion: Mitophagy plays an important role in the development of osteoporosis. The identification of three mitophagy-related genes may contribute to the early diagnosis, mechanism research and treatment of OP.

Keywords: bioinformatics; biomarker; differentially expressed genes (DEGs); mitophagy; osteoporosis; protein-protein interaction (PPI).

FIGURE 1

The flow chart of this research.

FIGURE 2

Detection of DEG between low and high BMD samples. (A) Volcano plot of differentially expressed genes identified in GSE56815 dataset. (B) Heat map of differentially expressed genes identified in GSE56815 dataset.

FIGURE 3

GSVA combined with WGCNA was used to screen MRG score related modules. (A) Sample clustering tree. (B) Data sample re-clustering and phenotypic information. (C) The soft threshold distribution. (D) The cluster dendregram after module merging. (E) The clustered modules.

FIGURE 4

Identification and functional enrichment of DE-MRGs. (A) The veen plot showed the interaction between DEGs and MRGs. (B) The top 10 functional enrichment in BP, CC, and MF analysis, respectively.(C) The KEGG of DE-MRGs.

FIGURE 5

OP-specific protein-protein interaction network.

FIGURE 6

Screening hub genes and evaluating their diagnostic value. (A) LASSO logic coefficient penalty diagram. (B) LASSO plot showed the variations in the size of coefficients for parameters shrank as the value of k penalty increased. (C) The relationship between the prediction accuracy of SVM-RFE and the number of features. (D) Boxplot of importance distribution of each gene in Boruta algorithm.

FIGURE 7

Evaluating the diagnostic value of hub genes (A) The veen plot showed the interaction of the LASSO, SVM-RFE and Boruta. (B) Receiver operating characteristic curve of three hub genes between OP group and control group. (C) nomogram model. (D) Calibration curve to evaluate the predictive ability of nomogram model.

FIGURE 8

The GESA of OP hub genes. (A) The GSEA of NELFB in OP. (B) The GSEA of SFSWAP in OP. (C) The GSEA of MAP3K3 in OP.

FIGURE 9

The IPA Pathway Analysis of OP hub genes. (A) Hub gene-related diseases and functional enrichment results. (B) MAP3K3 activation inhibition pathway diagram (left is the activation pathway, right is the inhibition pathway).

FIGURE 10

The immune cell infiltration association with hub genes. (A) Abundance box plots of 29 immune gene sets in OP group and Control group. (B) Heat map of the relationship between immune-related gene sets. (C) Bubble plot of correlation between OP hub genes and significantly different immune-related gene sets.**p < 0.01,***p < 0.001.

FIGURE 11

Network construction based on hub gene. (A) The veen plot of miRNA intersection predicted by miRwalk database and miRTarBase database. (B) The veen plot of lncRNA intersection predicted by ENCORI database and lncBase database. (C) CeRNA network diagram. (D) TF-mRNA network diagram.

FIGURE 12

Drug prediction and gene expression levels. (A) Interaction network diagram of hub genes and potential drugs. (B) Hub genes in the data set GSE56815 expression box plot.

FIGURE 13

RT-qPCR verification of hub genes. (A) NELFB. (B) SFSWAP. (C) MAP3K3. All results were expressed as mean ± standard deviation.(*p < 0.05; ns, not significant).