Yu Su1†
Gangying Yu2†
Dongchen Li1†
Yao Lu1Cheng Ren1Yibo Xu1Yanling Yang3Kun Zhang1
Teng Ma1*Zhong Li1*- 1Honghui Hospital, Xi’an Jiaotong University, Xi’an, China
- 2Department of International Ward (Orthopedic), Hospital of Chengdu University of Traditional Chinese Medicine, Chengdu, China
- 3Basic Medical College of Yan’an University, Yan’an, China
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.
1 Introduction
Osteoporosis (OP) is a systemic bone disease and the most common metabolic bone disease (Langdahl, 2020; Yan et al., 2023). They can be divided into two categories: primary and secondary. The former mainly occurs in children and young adults, considering the possible influence of genetic factors; however, its specific pathogenesis remains unknown. The latter is more common, mainly affecting older men and postmenopausal women, and is related to aging (Sozen et al., 2017; Guo et al., 2021). The main feature of OP is uncoupled bone resorption, which leads to decreased bone mass, damage, and degradation of the bone tissue microstructure, resulting in increased bone fragility and increased susceptibility to fractures, especially in the hip, wrist, and spine (Liu et al., 2021). OP is reported to be associated with fractures in more than half of white women and one-third of white men (Golob and Laya, 2015). Age is the most important risk factor of OP(6). In addition, estrogen deficiency, improper use of glucocorticoids, diabetes, malnutrition, and heavy drinking have been identified as important risk factors for OP(6, 7). The invisibility of bone loss and absence of symptoms at the beginning of OP make it a silent disease, as it is usually difficult to detect before the first fracture; however, damage can be fatal once it occurs (Unnanuntana et al., 2010; Johnston and Dagar, 2020).
Currently, the diagnostic methods for OP are based on Dual-Energy X-Ray Absorptiometry (DEXA) evaluation of bone mass (Chen et al., 2023), double X-ray absorption measurement (Aibar-Almazan et al., 2022), quantitative computed tomography, and bone tissue biopsy (Filippiadis et al., 2018). However, owing to its high detection cost, strong invasiveness, and large amounts of ionizing radiation, its long-term use is limited. There are three types of treatment options for OP, which are divided into anti-reabsorption therapy, anabolic therapy, and dual-acting therapy (Foessl et al., 2023). Anti-reabsorption therapy involves the use of anti-reabsorption drugs, including bisphosphonates, denosumab, selective estrogen receptor modulators (SERMs), and estrogen, to reduce bone resorption by inhibiting osteoclast function (Chen et al., 2023; Foessl et al., 2023). Anabolic drugs, such as teriparatide, which is used in anabolic therapy, is strong bone formation promoter effective in treating osteoporosis (Che et al., 2023). Romosozumab is currently the only drug available as dual-acting therapy. The dual effects of increased bone formation and reduced bone resorption lead to rapid increases in trabecular and cortical bone mass, which can improve bone structure and strength (Langdahl, 2020). However, most drugs have limitations and adverse side effects, and the current clinical treatments for OP are not ideal. Osteoclast-mediated bone destruction can occur, resulting in impaired therapeutic effects (Nogues and Martinez-Laguna, 2018). Therefore, there is an urgent need to develop new diagnostic methods for OP and provide new targets for the development of new drugs and immunotherapy by understanding the biological pathways of the markers.
Autophagy is a membrane-dependent, subcellular component turnover mechanism in eukaryotic cells. It is a complex dynamic process that maintains intracellular homeostasis under different physiological or pathological conditions by controlling protein degradation and the renewal of damaged organelles (Wang J. et al., 2023). Mammalian autophagy is divided into three subcategories: macroautophagy, microautophagy, and chaperone-mediated autophagy. Macroautophagy is the main autophagic process that helps effectively transport cytoplasmic cargo to lysosomes for degradation (Mizushima and Levine, 2020). Autophagy is highly regulated and an appropriate level of autophagy can protect cells from pathological or physiological damage. Conversely, excessive or reduced autophagy trigger apoptosis (Scheibye-Knudsen et al., 2014). Studies have shown that autophagy is necessary for cell survival. Dysregulation of autophagy is associated with various human diseases, including cancer and metabolic, neurodegenerative, cardiovascular, and age-related diseases (Mattson et al., 2008).
Mitochondrial autophagy is a selective form of autophagy. It is a autophagy process that occurs in mitochondria. This is a highly specific quality control process and the only known selective method for removing entire mitochondria. It maintains mitochondrial homeostasis by eliminating damaged organelles and redundant proteins, and reducing cell stress caused by harmful stimuli (Zong et al., 2016; Wang J. et al., 2023). Mitochondrial autophagy maintains cellular homeostasis under healthy conditions. However, it can also be induced under certain pathological or physiological conditions, thereby degrading healthy mitochondria and ultimately promoting the development of many diseases (Fang et al., 2019; Xu et al., 2020; Zhang et al., 2021). In addition, mitophagy has shown excellent medical application value in many chronic diseases and is considered to be an important target for the treatment of chronic diseases (Yan et al., 2023). Increasing evidence shows that abnormal mitochondrial autophagy can disrupt the balance of bone metabolism and plays a key role in bone metabolism disorders (Wang et al., 2017; Naik et al., 2019; Wang S. et al., 2020). Other studies have shown that mitophagy is involved in the regulation of OP. For example, phosphoinositide 3-kinase (PI3K) can treat or affect the progression of OP by interfering with mitophagy in bone marrow-derived mesenchymal stem/stromal cells (BMSC) (Yang et al., 2019; Chen et al., 2020). Lee found that the expression of Pten-induced putative kinase 1 (PINK1) was decreased in OP patients, and the PINK1/Parkinson‘s disease-related gene (Parkin) pathway is an important way of mitophagy (Lee et al., 2021). Maity et al. (2022) reported that the expression levels of mitophagy-related molecules PINK1 and Parkin in DPSC were increased, and mitophagy was induced to promote osteoblast differentiation. It has also been reported that epigallocatechin-3-gallate (EGCG) is a rich polyphenol in green tea, which can significantly reduce the expression of mitochondrial autophagy related PINK1 and Parkin. inhibit osteoclast differentiation, and regulae the progression of osteoporosis (Sarkar et al., 2022). All of the above prove that mitophagy may be related to OP. Therefore, we have reason to suspect that osteoporosis is associated with mitophagy. However, previous studies have only explored the effect of a specific gene-mediated mitophagy pathway on osteoporosis, suggesting that there is a certain correlation between the two. There is still a lack of research on the clinical significance of mitophagy-related gene (MRG) in osteoporosis. Therefore we conducted this study to understand the expression of mitophagy-related genes in OP, in order to provide an important reference for the diagnosis, mechanism research and treatment of OP.
2 Methods and materials
2.1 Data sources and downloaded data
The clinical data of all patients were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). According to the following criteria: (Yan et al., 2023) Disease name: Osteoporosis (Langdahl, 2020); Organism type: Homo sapiens (Guo et al., 2021). Sample source: blood sample (Guo et al., 2021). Dataset including samples with low and high bone mineral density (BMD) (Sozen et al., 2017). Sample size as large as possible (The dataset contains more than 3 samples per group). Hence, GSE56815, including the gene expression data (Expression profiling by array) of peripheral blood monocytes from 40 low-BMD samples and 40 high-BMD samples, was analyzed in this study. Thirty-four mitophagy-related genes (MRG) were downloaded from previous literature (Xu et al., 2022).
2.2 Identification of differentially expressed genes
We used the Wilcoxon test method in the R package to analyze the genetic data and get individual p-values after we completed the data download and collation (Cheng L. et al., 2022). The false discovery rate (FDR) was used to test the p-value, and the FDR value obtained was used to screen the differentially expressed genes (DEGs), with the following criterion:FDR<0.05 (Hsu and Lee, 2010). Volcano plots and heat maps were generated to identify the DEGs.
2.3 Screening of module genes related to mitophagy
In the GSE56815 dataset, the expression matrices of 34 MRG were extracted, and the MRG scores between OP and control samples were calculated using the gene set variation analysis (GSVA) package (version 1.42.0) (Suarez-Farinas et al., 2010). The weighted gene co-expression network analysis (WGCNA) package (version 4.0.3) (Langfelder and Horvath, 2008) in R was used to screen out module genes with high correlation with MRG score. Specifically, samples were clustered and outliers were excluded to ensure the accuracy of the analysis, and the appropriate soft threshold was determined with the following conventional criterion: 0.8 < R2 < 0.95, and the mean value of the adjacency function was close to 0. A co-expression matrix was constructed, and a hybrid dynamic clipping tree method with a minimum of 30 modules was used to identify co-expressed gene modules (Reyes et al., 2017), and the huge number of genes was classified into dozens of gene modules. Finally, the significant modules correlated with MRG scores were defined with a relevance |Cor|>0.3, p < 0.05.
2.4 Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analysis of differentially expressed-MRG
The online jvenn website (http://jvenn.toulouse.sra.fr/app/example.html) was used to intersect the module MRG screened by DEGs and WGCNA to obtain differentially expressed-MRG, which were denoted as DE-MRG. The R package “clusterProfiler” (version 4.4.4) (Yu et al., 2012) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DE-MRG. A p-value <0.05 was considered statistically significant. The R package “ggplot2” (version 3.3.2) (Ito and Murphy, 2013) was used to visualize the enrichment results.
As part of the GO analysis, three ontologies were identified, namely molecular function (MF), cellular component (CC), and biological process (BP), which helped to explore the biological processes of these DEGs. KEGG analysis was used to identify metabolic or signal transduction pathways that were significantly enriched in the DEGs.
2.5 Specific protein-protein interaction network construction
At present, protein-protein interaction (PPI) networks are among the best-studied biomolecular networks. To explore the interaction between DE-MRG, we used the STRING (https://string-db.org) website to construct the PPI network of DE-MRG and screen the threshold: confidence = 0.4 (Zhao et al., 2015), and then visualized with Cytoscape (Shannon et al., 2003).
2.6 Signature gene identification
We identified candidate key genes by the intersection of key module genes and DGEs. Subsequently, three machine learning algorithms, least absolute shrinkage and selection operator (Lasso), support vector machine recursive feature elimination (SVM-RFE), and Boruta were used to analyze the hub genes.
Lasso regression is a new variable selection technique, which is implemented by using the glmnet package with penalty parameters, and the feature genes are obtained after 10-fold cross validation. SVM-RFE is an effective feature selection technique using the “e1071” package and “glmnet” package (Friedman et al., 2010) to obtain the relationship between the prediction accuracy and the number of features, and the number of hub genes is obtained when the generalization error is the lowest. The Boruta algorithm is a wrapper method built around a random forest classifier (Shen et al., 2021). This algorithm minimizes the error of the random forest model, which eventually forms a subset of the minimum optimal features.
Using the online Jveen map-making website (http://jvenn.toulouse.inra.fr/app/example.html), overlapping genes were obtained by intersecting the markers derived from the three algorithms, that is, the key genes obtained by screening.
The receiver operating characteristic (ROC) curve of the hub genes was drawn using the R pROC package (version 1.12.1) (Hanley and McNeil, 1982) and the diagnostic efficiency of the key genes was evaluated using the area under the curve (AUC). An AUC greater than 0.7 indicated good diagnostic performance.
The nomogram model can evaluate the relationship between variables in the prediction model (Park, 2018). The nomogram model of the hub genes was constructed using the R RMS package (version 6.0-1) (Zhao and Li, 2023), and a calibration curve was drawn using the calibration function and boot method. The accuracy of the nomogram model for predicting osteoporosis was determined using the slope of the curve and Hosmer-Lemeshow goodness-of-fit test (HL test), where the insignificant difference (p > 0.05) between predicted and actual observed results in HL test indicates that the prediction model has good calibration ability (Wang et al., 2021).
2.7 Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analysis of hub genes
The Pearson correlation coefficient between the key genes and the remaining genes was calculated and ranked from highest to lowest. After the corresponding gene sequence list of each key gene was obtained, the “clusterProfiler” package (version 4.4.4) (Yu et al., 2012) was used for GO and KEGG Gene Set Enrichment Analysis (GSEA), and Benjamini & Hochberg method was used for multiple test correction. The corrected p-value was p.adjust, and p.adjust <0.05 was considered as significant enrichment.
2.8 Classical pathways of the DEGs were analyzed by IPA pathway analysis
IPA Pathway Analysis (IPA) is an integrated analysis software. We used DEGs to perform IPA functional enrichment analysis to identify typical pathways, diseases, and functional pathways related to key genes. A z-score >2 indicates that the pathway is activated, and a z-score < -2 indicates that the pathway is inhibited. All pathways containing key OP genes were screened for display.
2.9 Immune microenvironment analysis
Single sample GSEA (ssGSEA) was performed by calculating the rank value of each gene based on the expression profile. To explore the composition of immune cells in the GSE56815 dataset, ssGSEA package (version 1.36.3) (Zuo et al., 2020) was used to perform immune cell infiltration analysis, and the differentially expressed immune cells between the disease group and the control group were detected using Wilcoxon test (p value less than 0.05 was considered considerable). The results were presented as boxplots. Subsequently, the “spearman” algorithm was used to calculate the correlation between the immune-related gene sets and the correlation between OP hub genes and the significantly different immune-related gene sets, respectively. The results were presented as heat-and-bubble maps.
2.10 Network construction based on hub gene
2.10.1 Construction of the ceRNA network
For the obtained OP hub genes, we used the miRWalk database (http://mirwalk.umm.uni-heidelberg.de/search_genes/) and miRTarBase database (https://mirtarbase.cuhk.edu.cn/∼miRTarBase/miRTarBase_2022/php/search.php) to search for relevant miRNAs, then took the intersection to obtain the miRNA of OP hub genes. Based on the obtained miRNAs, the databases ENCORI (https://starbase.sysu.edu.cn/agoClipRNA.php?source=mRNA) (the default parameters are set: CLIP-Data ≥ 1, Degradome-Data ≥ 0, pan-Cancer ≥ 0) and LncBase (https://diana.e-ce.uth.gr/lncbasev3/interactions) were searched; then, the corresponding lncRNAs of the miRNA obtained in the two databases were intersected to obtain the lncRNA of OP hub genes. Cytoscape (Shannon et al., 2003) was used to draw the visualized mRNA-miRNA-lncRNA network diagram.
2.10.2 Construction of TF-hub gene network and predicting potential drugs
Based on the obtained key OP genes, NetworkAnalyst (https://www.networkanalyst.ca/NetworkAnalyst/uploads/ListUploadView.xhtml) was used to search for the TF corresponding to the key genes to obtain the interaction information between them, and parameter settings are as follows: Specify organism = “H. sapiens (human)”, Set ID type = “Official Gene Symbol” Cytoscape (Shannon et al., 2003) was used to visualize the interactions.
2.11 Predicting potential drugs
Potential drugs for the treatment of OP were predicted using the CTD database (http://ctdbase.org/), and the interactions between the two were analyzed. Cytoscape (Shannon et al., 2003) was used for the visual display.
2.12 Validation of key gene expression
From the GSE56815 dataset, the expression levels of key OP genes were extracted and analyzed using the Wilcoxon test. The expression of key genes was checked and the results are displayed as box plots.
2.13 Experimental validation by RNA extraction and real-time quantitative polymerase chain reaction
The study was in accordance with the Declaration of Helsinki (as revised in 2013) (World Medical Association, 2013). This study was approved by the Ethics Committee of Honghui Hospital Affiliated to Medical College of Xi ‘an Jiaotong University (No.202212012). All participants were informed and signed consent. The study included patients who met the criteria in the orthopedics department of the hospital from January 2023 to November 2023.
Inclusion criteria: (Yan et al., 2023) Those who met the diagnostic criteria of OP, according to “Chinese expert consensus on the diagnosis of osteoporosis by imaging and bone mineral density” (Cheng et al., 2020; Langdahl, 2020). All subjects were informed and voluntarily signed the consent form.
Exclusion criteria (Yan et al., 2023): those who do not meet the diagnosis of OP; (Langdahl, 2020) patients with severe cardiovascular disease (Guo et al., 2021); Patients with severe liver/kidney dysfunction (Sozen et al., 2017); patients with coagulation dysfunction (Liu et al., 2021); patients with autoimmune diseases (Golob and Laya, 2015); untreated bone metabolic diseases other than osteoporosis (Johnston and Dagar, 2020); Long-term use of drugs affecting bone metabolism (such as glucocorticoids and steroid hormones) (Unnanuntana et al., 2010); Anti-osteoporosis treatment history (such as bisphosphonates, teriparatide and calcitonin, etc.).
Grouping criteria: lumbar spine and total hip BMD (DEXA): normal group T-score≥1.0; T-score ≤ -2.5 in OP group (all BMD measurements were performed by the same trained professional).
Finally, 8 patients with osteoporosis and 8 patients with normal bone mineral density were included in the study. There were 4 males and 4 females in each group. Peripheral venous blood (10 mL) was collected from all patients after fasting for one night (≥12 h).
The operation of human peripheral blood mononuclear cells (PBMC) is carried out as described above (Chen et al., 2021). First, the patient‘s whole blood was placed in a 50 mL centrifuge tube, diluted with 10 mL PBS and gently mixed. Subsequently, the blood sample was continuously centrifuged for 20 min at a speed of 2000 rpm. After centrifugation, the white blood cell layer in the center of the sample containing PBMC was sucked out with a pipette and transferred to a new 15 mL centrifuge tube, and 10 mL PBS was added to it. The centrifuge was used to continuously centrifuge for 10 min at a speed of 1500 rpm, and then the supernatant was taken out for precipitation, which was the required PBMC. The obtained PBMC were inoculated in a 6-well plate, and 1 mL of TRIzol (GENSTAR Inc. Beijing, China) reagent was added to each well to extract total RNA from the cells. Subsequently, 1 mg total RNA was reverse transcribed using a cDNA synthesis kit (Thermo Fisher Scientific Inc. Shanghai, China). The cDNA was detected using a 20 mL SYBR Green qPCR SuperMix (TargetMol Chemicals Inc. Shanghai, China) and RT-qPCR machine (Bio-Rad, Hercules, CA, United States). The final thermal cycle conditions for gene amplification were: 95 C 30 s, 95 C 5 s 40 cycles, and the last step was 60 C 30 s. Quantitative analysis was performed using the 2△△CT method to calculate the relative expression of each gene. The hub gene-related detection primers were prepared by Beijing Tsingke Biotech Co., Ltd. (China), as shown in Table 1.
TABLE 1. Sequences of hub gene-specific primers used for RT-qPCR.
2.14 Statistical analysis
All statistical analyses were performed using R software. Differentially expressed analysis for DEGs, differentially expressed immune cells, and key OP genes between OP and control samples in GSE56815 datasets were conducted using Wilcoxon test. When p < 0.05, all results were considered statistically significant and all p-values were two-tailed.
The flow chart of this study is shown in Figure 1.
FIGURE 1. The flow chart of this research.
3 Results
3.1 Differential gene analysis
In the GSE56815 dataset, 548 genes were detected between low- and high-BMD samples, of which 366 genes were upregulated and 182 genes were downregulated (Figures 2A, B; Supplementary Table S1).
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.
3.2 GSVA and WGCNA screen mitochondrial autophagy-related genes scoring related modules
The expression matrix of the 34 MRG was obtained from the dataset GSE56815. MRG scores were calculated between OP samples and control samples using the “GSVA” package.
The WGCNA package was used to analyze GSE56815, and the MRG score was used as WGCNA trait data to screen out module genes with high correlation with the MRG score; the samples were clustered to show the overall correlation of the samples in the data set. The threshold was set to 42 and outlier samples were eliminated, that is, those above the red line in the figure (Figure 3A). It was then reclustered, and the corresponding sample traits were added (Figure 3B). When β = 11, the co-expression network most closely approximates the scale-free distribution (Figure 3C). When MEDissThres = 0.3 was set, similar modules were analyzed, and the number of modules that aggregated after merging was 16 (Figure 3D). Analysis of the correlation of each module and the MRG score of red represents MRG scores are related, green represents is negatively related to the MRG score, color represents the degree of correlation, in | Cor | > 0.3, p < 0.05, for threshold selection, obtain results with MRG score high correlation modules with a total of three, blue (number of genes:1947), pink (number of genes:278), grey60 (number of genes:66), and the total number of genes was 2291(Figure 3E).
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.
3.3 Identification and functional enrichment of differentially-expressed mitochondrial autophagy-related genes
The 2291 hub module genes obtained by WGCNA overlapped with 548 DEGs obtained by differential analysis between the high- and low-BMD groups, and a total of 91 overlapping genes were obtained, namely, DE-MRG (Figure 4A). DE-MRG were subjected to GO functional enrichment and KEGG pathway analysis and were enriched in 430 biological processes (BPs), 59 molecular functions (MFs), 40 cellular components (CCs) (Figure 4B), and 33 pathways (Figure 4C) (p < 0.05). For BPs, the top 10 enrichment programs were mainly involved in the positive regulation of response to external stimuli, I-kappaB kinase/NF-kappaB signaling, cellular response to biotic stimulus, homeostasis of number of cells, import into nucleus, and other biological processes. For CC and MF, the top 10 enriched items are shown in Figures, respectively, and the enrichment analysis of KEGG pathway showed that DE-MRG are primarily involved in Salmonella infection, Fc gamma R-mediated phagocytosis, Epstein-Barr virus infection, and neurotrophin signaling pathways.
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.
3.4 Protein-protein interaction analysis
The construction of PPI network showed that gene AKT1 was the central gene of this PPI network, and it interacted with the most proteins. The results were visualized using the Cytoscape software (Figure 5).
FIGURE 5. OP-specific protein-protein interaction network.
3.5 Screening for the best characteristic genes of osteoporosis and evaluating their diagnostic value
A total of 18, 7, and 28 genetic biomarkers were identified using the machine learning algorithms LASSO, SVM-RFE, and Boruta, respectively (Figures 6A–D). Moreover, the key genes shared by the three OP were identified by overlapping biomarkers derived from the three algorithms, including negative elongation factor complex member B (NELFB), Splicing Factor (SFSWAP), and mitogen-activated protein kinase kinase 3 (MAP3K3) (Figure 7A).
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.
We determined the reliability of the diagnostic value of our hub genes by plotting a ROC curve and calculating the area under the AUC; those of NELFB, SFSWAP, and MAP3K3 were 0.75(95% CI = 0.64-0.86), 0.71(95% CI = 0.59-0.82), and 0.70(95% CI = 0.58-0.81), respectively, with all greater than 0.7 (Figure 7B). Therefore, our gene hub can effectively diagnose osteoporosis. After the construction of the nomogram model of the gene hub(C-index = 0.886, 95% CI = 0.815-0.957), the slope of the calibration curve was close to 1, and the HL test exhibited satisfactory agreement between predicted and actual outcomes with p = 0.545 > 0.05, indicating that the nomogram model had high prediction accuracy for osteoporosis (Figures 7C, D).
3.6 Functional enrichment analysis of hub genes
GSEA enrichment analysis of GO and KEGG was performed on the three hub genes, and p.adjust <0.05 was considered a significant enrichment result (Supplementary Table S2). According to the significance ranking, GO and KEGG TOP5 Descriptions were selected for visual display (Figures 8A–C).
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.
We evaluated the signaling pathways associated with the key genes using GSEA. The first five signaling pathways are shown in the figure. GO analysis showed that NELFB was significantly associated with cytoplasmic translation, cytosolic, large ribosomal subunit, cytosolic ribosome, ribosome, and structural constituents of the ribosomes. SFSWAP expression was significantly correlated with cytoplasmic translation, cytosolic ribosomes, ribosomal subunits, ribosomes, and the structural constituents of ribosomes. MAP3K3 expression significantly correlated with cytoplasmic translation, cytosolic ribosomes, ribosomal subunits, ribosomes, and structural constituents of ribosomes.
KEGG results showed that NELFB was significantly correlated with the B cell receptor signaling pathway, chemokine signaling pathway, Fc gamma R-mediated phagocytosis, human cytomegalovirus infection, and ribosomes. Sfswap expression was significantly associated with Fc gamma R-mediated phagocytosis, oxidative phosphorylation, Parkinson’s disease, ribosomes, and spliceosomes. The expression of MAP3K3 was significantly associated with Fc gamma R-mediated phagocytosis, osteoclast differentiation, oxidative phosphorylation, Parkinson’s disease, and ribosomes.
In summary, GO analysis showed that the molecular functions of these genes were mainly enriched with cytoplasmic translation, cytosolic ribosome, ribosome and structural constituent of ribosome; KEGG pathway analysis showed that it was mainly enriched in Fc gamma R-mediated phagocytosis and Ribosome.
3.7 IPA pathway analysis
The results of the IPA functional enrichment analysis showed 20 inhibition pathways (blue) and activation pathways (range) (Figure 9A). The hub gene, MAP3K3, was included in these 40 pathways. The genes were involved in six pathways, and three activation pathways were GNRH Signaling, Natural Killer Cell Signaling and G-Protein Coupled Receptor Signaling. The three inhibited pathways were Cardiac Hypertrophy Signaling (enhanced), RANK Signaling in Osteoclasts, and Cardiac Hypertrophy Signaling. The strongest activation intensity of MAP3K3 (GNRH Signaling) and strongest inhibition intensity of the pathway (cardiac hypertrophy Signaling) are shown in Figure 9B.
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).
3.8 Immune infiltration analysis of the dataset
Immune infiltration analysis was performed on the integrated GSE56815dataset. The scores of the 29 immune cell types evaluated using ssGSEA showed significant differences in the expression of four immune-related gene sets: APC_co_stimulation, iDCs, Tfh, and Th1_cells (Figures 10A, B) (Supplementary Table S3). The correlations results between immune-related gene sets are shown in Figure 10B. (Supplementary Table S4). Further, the heatmap for the correlation of hub genes and immune-related gene sets suggested that there was a significant negative correlation between NELFB and iDCs in hub genes with significant differences (cor = −0.390, p < 0.001). There was a positive correlation with Tfh levels (cor = 0.233, p < 0.05) (Figure 10C).
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.
3.9 Network construction based on hub genes
We identified eight miRNAs, of which seven corresponded to MAP3K3 and one to NELFB (Figure 11A). Simultaneously, 81 lncRNAs were obtained, corresponding to four miRNAs (hsa-miR-132-3phsa-miR-182-5p, hsa-miR-212-3p, and hsa-miR-324-5p), which correspond to the key genes MAP3K3 and NELFB (Figure 11B). Based on the above two intersection pairs, 79 pairs of lncRNA-miRNA-mRNA relationships regulated by the same miRNA were screened, including 68 lncRNAs, four miRNAs, and two mRNAs, and a ceRNA network was constructed, that is, ceRNA network (Figure 11C).
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.
In addition, 72 TF–mRNA interaction pairs were retrieved, and only the key genes, SFSWAP and MAP3K3, provided TF interaction information. Visual expression was performed using Cytoscape software (Figure 11D).
3.10 Potential drugs prediction
By predicting potential drugs for the treatment of OP in the CTD database, 17 drug-gene pairs were obtained. Both vinclozolin and bisphenol A are associated with SFSWAP and MAP3K3. Benzo (a) pyrene and valproic acid were associated with both SFSWAP and NELFB. The results were visualized using Cytoscape software (Figure 12A).
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.
3.11 Dataset validation
The box plot shows that in the dataset GSE56815, compared with the control group, the genes MAP3K3 and SFSWAP were significantly expressed (p < 0.01), whereas the expression of NELFB was lower (Figure 12B).
3.12 RT-qPCR verification of hub genes
The results of RT-qPCR showed that the mRNA expression level of NELFB in OP group was lower than that in normal group. There was no difference in the mRNA expression levels of SFSWAP and MAP3K3 between the OP group and the normal group. The expression level of hub gene NELFB was consistent with bioinformatics (Figures 13A–C).
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).
Generally speaking, Wilcoxon test combined with WGCNA methods were used to screen 91 DE-MRG related to the MRG score, in which AKT1 interacted with the most proteins in the PPI network. Next, NELFB, SFSWAP, and MAP3K3 were selected as the key genes with excellent diagnostic value for OP. These key genes was mainly enriched in the biological terms of cytoplasmic translation and Fc gamma R-mediated phagocytosis, especially, MAP3K3 involved in the activation and inhibition of various IPA Pathway, such as Cardiac Hypertrophy Signaling and GNRH Signaling. Moreover, there were considerable correlations between NELFB and iDCs (cor = −0.39), and Tfh levels (cor = 0.23). It is worth noting that the expression level of NELFB detected by RT-qPCR was consistent with the results of bioinformatics. Finally, the potential ceRNA network, TF binding sites, drugs targeting three key genes was predicted. For the clinical utilize, the nomogram constructed by converting the expression of three hub genes into a total score had a high prediction accuracy for osteoporosis (C-index = 0.886, 95% CI = 0.815-0.957; HL test p > 0.05).
4 Discussion
Basic research on mitophagy and osteoporosis has received considerable attention. This study evaluated the DEG between patients with OP and healthy cohorts and explored the key modules of MRG based on GSVA-WGCNA. The intersection of the 91 DE-MRG was analyzed using KEGG and GO. A PPI network was constructed, and the core target gene was AKT1. Three hub genes related to OP and mitophagy genes were identified by three machine learning methods, including NELFB, SFSWAP, and MAP3K3. ROC and nomogram models proved their diagnostic value. GO and KEGG analyses showed that the three hub genes were involved in the ribosome and cytoplasmic ribosome pathways. IPA showed that MAP3K3 was associated with six pathways, including GNRH Signaling. The ssGESA indicated that NELFB was highly correlated with iDCs (cor = −0.390, p < 0.001). The network based on hub genes showed that 3 TFs (TRIM22, CEBPG, and CEBPZ), 4 miRNAs (hsa-miR-132-3p, hsa-miR-182-5p, hsa-miR-212-3p, and hsa-miR-324-5p), and 4 drugs were involved in osteoporosis. Finally, the expression level of NELFB detected by RT-qPCR was basically consistent with bioinformatics.
We performed DE-MRG enrichment analysis on the selected samples, and the results showed that they were involved in many biological processes, including “I−kappaB kinase/NF−kappaB signal”, “regulation of podosome assembly”, and “podosome assembly”. The inhibition of NF−kappaB signal can affect the formation of osteoclasts and has a high correlation with osteoporosis. In addition, the most common molecular functions include ‘heterochromatin’ and ‘secretory granule lumen’. Aging of BMSC is the main cause of osteoporosis. Several studies have reported that related genes can promote or inhibit the occurrence of osteoporosis by regulating heterochromatin to interfere with MSC(51). For example, KDM4B inhibits the formation of H3K9me3 (heterochromatin) lesions to reduce MSC self-renewal and affect osteoporosis (Deng et al., 2021), whereas sIRT3 inhibits the generation of heterochromatin through the vulcanization of hydrogen sulfide, thereby inhibiting the senescence of BMSC and providing a new target for the treatment of osteoporosis (Liu et al., 2023). Secretory granules indirectly affect blood glucose levels by regulating insulin secretion, and diabetes is a clinical risk factor for OP.
The central node of the PPI molecular network is AKT1, which interacts with and performs molecular docking between small molecules and encoded proteins. AKT1 is a serine/threonine protein kinase known as Akt kinase (Akt1, Akt2 and Akt3) (Xia et al., 2020). Previous studies have shown that induced AKT1 expression promotes the proliferation of mesenchymal stem cells and ultimately inhibits their apoptosis, thereby alleviating osteoporosis (Wang et al., 2022). In addition, AKT1, as a core target gene, can regulate OP through the PI3K-Akt signaling pathway, which is manifested in the phosphorylation of AKT1 and the expression of PI3K to promote osteogenesis and regulate the progression of OP (Yang et al., 2019; Zhao et al., 2021; Liang et al., 2022). AKT has also been reported to affect osteoporosis by upregulating FOXO1 and enhancing the expression of bone turnover markers (ALP, OCN, Runx2, and Col1) and extracellular matrix mineralization (Cai et al., 2021). AKT1 plays a key role in OP by participating in multiple signaling pathways.
NELFB is a negative elongation factor B, also known as COBRA1, which is one of the four subunits of the NELF complex that controls the participation of RNA polymerase II (Pol II) in gene transcription and is essential for multiple biological processes. Osteoporosis is caused by an imbalance between osteoclasts and osteoblasts and occurs near hematopoietic cells in the bone marrow. Osteoblasts are the primary cells involved in bone formation and play a central role in hematopoiesis (Salamanna et al., 2020; Kim et al., 2021). Researchers have found that the absence of NELFB can induce excessive progenitor cell development during primitive and final hematopoiesis. Hematopoietic progenitor cells can bind hematopoietic growth factors to promote interactions between osteoblasts and hematopoietic stem cells (HSC), thereby affecting the number and activity of osteoblasts (Huang et al., 2022). Therefore, NELFB and osteoblasts appear to be associated with osteoporosis via hematopoietic progenitor cells (Calvi et al., 2003; Kim et al., 2011). It has been reported that NELFB has a synergistic effect with TCF1 in T cell response to cancer, while a non-coding RNA AK045490 inhibits osteogenic differentiation by down-regulating the expression of TCF1, LEF1 and Runx2. NELFB may affect osteogenic differentiation and bone formation by interfering with TCF1, thus providing a potential new drug target for osteoporosis (Li et al., 2019). In addition, NELFB can inhibit the transcription of estrogen receptor (ER) receptor gene, interfere with its binding to estrogen, inhibit the expression of osteoprotegerin (OPG), and promote the role of nuclear factor (NF)-κB ligand (RANKL), thereby promoting osteoclast formation and bone resorption activity, which is closely related to the pathogenesis of osteoporosis (Cheng CH. et al., 2022).
SFSWAP is a hypothetical splicing factor that encodes a protein containing an RS domain (SR-like). Proteins containing RS domains regulate RNA processing, including splicing, transcript extension, transcript stability, nuclear export, miRNA cleavage, and genomic stability (Moayedi et al., 2014). By studying the mutation process of the Notch signaling pathway, Moayedi et al. found that the SfSWAP gene has a synergistic relationship with the sub-allele of the Notch ligand Jagged1 (Jag1), which can play the same role in the same genetic pathway and interfere with other genes of Notch signaling (Moayedi et al., 2014). We speculate that SfSWAP is associated with Notch signaling pathway and Jag1 gene.
JAG1 is a traditional osteoporosis gene, which was confirmed using the cFDR method in a study by Hu et al. (2018). Through a large genome-wide association study and a follow-up replication study, Gong et al. identified JAG1 as a candidate gene for BMD regulation in different races and as a potential key factor in the pathogenesis of fractures (Kung et al., 2010). Researchers signaling has been reported to strictly control mammalian cell fate determination during embryogenesis and adulthood, which is essential for skeletal development and activity of skeletal cells (Zanotti and Canalis, 2013). Disorders in Notch signaling are related to human diseases that affect bones. Wang et al. identified JAG1 and NOTCH1 as L-R genes with ossification-related functions using single-cell RNA sequencing, which further confirmed the inevitable relationship between the two and osteoporosis (Wang S. et al., 2023). Through animal experiments, Hui et al. found that the targeted inhibition of JAG1/NOTCH1 signal transduction can promote poor osteogenic differentiation of bone marrow mesenchymal stem cells (Han et al., 2022). These studies reflect the inevitable connection between these two factors and osteoporosis. In addition, considering that both Notch pathway and Jagged1 gene are associated with SfSWAP, this may indicate the potential role of SFSWAP in the development of osteoporosis.
MAP3K3, also known as mitogen-activated protein kinase kinase kinase kinase kinase kinase 3, is a member of the serine/threonine protein kinase family that is generally expressed and acts as an oncogene (Zhang Y. et al., 2019). Many studies have reported that MAP3K3 is highly correlated with the progression and treatment of nasopharyngeal carcinoma (Yin et al., 2019), ovarian cancer (Zhang Y. et al., 2019), lung cancer (He et al., 2015), and breast cancer (Fan et al., 2014). This is a universally expressed cancer-related gene. Another study has reported that MAP3K3 may be associated with cerebral cavernous malformations. It are believed that the mutation of MAP3K3 can define a subclass of cerebral cavernous malformations (Weng et al., 2021). In addition, He et al. found that six genes, including MAP3K3, may play an important role in the regulation of BMD in women by analyzing the association between multiple gene expression profiles and bone mineral density variation. This study also provides new therapeutic targets for the treatment of osteoporosis (He et al., 2016). However, there is no direct evidence of a link between MAP3K3 and OP. Perhaps it can be used as an object to study the molecular mechanism of OP in the future.
In addition, GSEA enrichment of the hub genes revealed that all three genes were involved in the ribosome and cytoplasmic ribosome pathways. Studies have found that endoplasmic reticulum stress is significantly correlated with osteoblast and osteoclast differentiation and osteoclast formation during osteoporosis progression. During this process, the number of ribosomes in the endoplasmic reticulum membrane decreases significantly (Li et al., 2017). This indicates that ribosomes are associated with the progression of osteoporosis. In addition, it has been reported that the pathogenesis of osteoporosis in the elderly may be related to ribosome-related genes and pathways (Wang C. et al., 2023). These results suggest that ribosomes and cytoplasmic ribosomes play key roles in osteoporosis.
The use of multiple drugs, including gonadotropin-releasing hormone (GnRH) agonists and aromatase inhibitors, increases the risk of fractures (Krishnan and Muthusami, 2017). Studies have reported that long-term treatment with gonadotropin-releasing hormone (GnRH) agonists induces competitive and reversible GnRH receptor blockade, thereby inhibiting the release of gonadotropins and sex hormones. A reduction in sex hormones can lead to a variety of adverse reactions, including accelerated bone loss, which positively promotes the occurrence of OP (Mohamad et al., 2021). This suggests that GnRH Signaling may be associated with osteoporosis. Natural Killer Cell Signaling has also been speculated to be related to OP. Osteoporosis is the result of a ‘bone remodeling’ imbalance, and there is a close relationship between the immune system, bone physiology, and pathology. Srivastava et al. developed the term ‘immune osteoporosis’ to emphasize the role of immune cells in osteoporosis pathology (Srivastava et al., 2018). Saxena et al. have reported that natural killer cells belonging to the lymphoid lineage have congenital characteristics and play a role in osteoporosis. There are also literatures that have concluded through bioinformatics methods that there are significant differences in natural killer cells between osteoporosis patients and non-OP patients (Saxena et al., 2021). This is consistent with our findings. GPR48, a member of the G protein-coupled receptor (GPCR) superfamily, has been reported to inhibit osteoclast differentiation by antagonizing the interaction between RANK and its ligand-RANKL, thereby interfering with the OP process (Filipowska et al., 2022). In addition, GPR125 positively regulates osteoclasts through RANKL-stimulated MAPK and AKT-NF-κB signaling pathways to participate in the treatment of osteoporosis (Tang et al., 2022). In conclusion, G-Protein Coupled receptor Signaling is highly correlated with the development, diagnosis, and treatment of OP. There are few reports on the relationship between the other three inhibited pathways and OP, and further research is required.
iDCs showed the highest correlation with key genes among the four immune-related genes based on ssGSEA. Studies have shown that dendritic cells (DCs) are effective antigen-presenting cells that are widely distributed in the bone immune and/or mucosal mesenchymal interface and can affect OP by activating RANKL-nuclear factor kB receptor activator or RANK-osteoprotegerin (Liu and Teng, 2023). In addition, we found a positive correlation between NELFB and Tfh levels (p < 0.05). However, there are few reports in this field, and further studies are required to confirm this.
Both mRNA and miRNAs are indispensable components of ceRNA regulatory networks. It is undeniable that the associated miRNAs are equally important (Xu et al., 2016). Identifying the relationship between miRNAs and diseases not only improves the understanding of molecular mechanisms and disease pathogenesis, but also facilitates clinical diagnosis and treatment (Zhang L. et al., 2019). In our study, four miRNAs, namely hsa-miR-132-3p, hsa-miR-182-5p, hsa-miR-212-3p, and hsa-miR-324-5p, were identified in the circRNA-miRNA-TF mRNA regulatory network. To date, no study has directly explored the roles of these four miRNAs in osteoporosis. However, our data further support the possibility that these four miRNAs are important regulators of osteoporosis development. Nevertheless, this conjecture and its precise mechanism need to be further studied.
A transcription factor (TF), also known as gene promoter, is a protein that can bind to gene-specific sequences to mediate gene transcription and expression (Sebastian and Contreras-Moreira, 2013). Various transcription factors play regulatory roles in OP pathogenesis (Liu et al., 2019). In this study, TF-mRNA relationship pairs, including SFSWAP and MAP3K3, were obtained. The transcription factor TRIM22 is related to both genes. Many studies have shown that the absence of TRIM22 can affect the PI3K/Akt/mTOR pathway, and studies have reported that the latter has a potential relationship with OP progression (Li et al., 2018). In addition, Jiang et al. reported that the expression of the CCAAT/enhancer-binding protein (CEBP) homologous protein can promote BMP4-induced osteogenesis of MSCs in vitro and in vivo. Transcription factors CEBPG and CEBPZ may have the same effect (Jiang et al., 2018). The TF-mRNA relationship that we studied may provide strong evidence for verifying the potential relationship between hub genes and OP.
CTD is a powerful public database and its significance lies in increasing our understanding of how environmental exposure affects human health (Davis et al., 2021). We used key genes combined with the CTD database to predict potential drugs for OP treatment. The results showed that vinclozolin and bisphenol A were simultaneously associated with SFSWAP and MAP3K3, and benzo (a) pyrene and Valproic Acid were associated with both SFSWAP and NELFB. Studies have reported that some drugs may lead to osteoporosis. For example, bisphenol A inhibits osteoblast differentiation and bone formation by activating RORα, leading to the formation of osteoporosis (Maduranga et al., 2022). Previous studies have shown that bisphenol A is associated with a high prevalence of lumbar osteopenia and osteoporosis in postmenopausal women (Wang N. et al., 2020). Valproic Acid activates the Notch signaling pathway and has a positive effect on bone defect repair. A previous study reported that the Notch signaling pathway is closely related to osteoporosis (Tao et al., 2021). There is no evidence in previous reports that other drugs, such as vinclozolin and benzo (a) pyrene, play a role in OP; therefore, our study first speculated that vinclozolin and benzo (a) pyrene may have a therapeutic effect on OP, which needs to be verified by further related experiments.
We note that two recent articles are similar to our research, and we analyze the similarities and differences between them. These articles are roughly the same as our research in the overall research method, and the subjects of the study are osteoporosis. But the biggest difference is the choice of genes. Song Hao et al.selected immune-related genes. Because bone and immune system have a common developmental niche, the development and remodeling of bone are affected by the immune system. Immune cells can participate in the pathogenesis of OP by producing pro-inflammatory mediators. Huang Xinzhou directly selected differentially expressed genes in osteoporosis as the research object. We selected mitophagy-related genes for analysis. Secondly, both articles use machine learning methods to identify genes. In the former study, LASSO and mSVM-RFE methods are selected. The latter was analyzed by LASSO and the Gaussian mixture model. In this study, three machine methods, Lasso, SVM-RFE and Boruta, were selected for screening. Therefore, although these two studies are similar to our research, the final screening genes and results are not the same. In summary, our study does not repeat the previous steps, and there is a lot of innovation.
The data in this study were all from the selected data set GSE56815, which was only produced by human peripheral blood monocytes (PBMC). A number of studies have found that PBMC are involved in the occurrence of OP (Liu et al., 2015; Xie et al., 2023). It can produce a variety of cytokines and growth factors that affect bone metabolism, such as macrophage colony stimulating factor, interleukin 1, IL-6 and transforming growth factor β (Fischer and Haffner-Luntzer, 2022; Xie et al., 2023). In addition, monocytes are also precursors of osteoclasts with bone resorption activity. Under certain conditions, osteoclasts are produced by monocytes in vitro (Sprangers et al., 2016). However, although monocytes are the premise of osteoclasts, bone marrow can include osteoblasts, osteoclasts and osteocytes. Studies have found that almost all pathways in PBMC overlap with pathways in bone marrow tissue (Xie et al., 2023). Therefore, peripheral blood mononuclear cells can only reflect the situation in the bone marrow to a certain extent, and the study of OP through PBMC is limited. This is one of the limitations of this study. In addition, other limitations that should be recognized include: the sample size of data from public databases is relatively small, which can easily lead to selection bias; the verification step only selected clinical trial verification, and the corresponding database was not selected for further verification, which reduced the credibility of the results.
Although several mitophagy-related genes were screened in this study, and the reliability of our analysis was verified from the side by biological analysis, several conclusions were also speculated, including that SfSWAP affected the progress of OP by association with Notch signaling pathway and Jag1 gene. The 4 miRNAs may be important regulators in the development of osteoporosis. Vinclozolin and Benzo (a) pyrene may have therapeutic effects on OP. However, the potential regulatory mechanism of mitophagy in the progression of OP has not been fully elucidated, and the conclusions we speculate have not been determined. In the future, more in vivo and in vitro experiments and biological analysis are needed to further verify.
5 Conclusion
In general, this study studied the association between mitophagy-related genes and OP through a number of machine learning and biological analysis, and screened three hub genes, NELFB, SFSWAP and MAP3K3, which can be used as potential biomarkers for the diagnosis, mechanism research and treatment of osteoporosis. In addition, our study supports a more in-depth study of the mechanism of action of these three mitophagy-related genes in OP.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Ethics statement
The studies involving humans were approved by the Medical Ethics Committee of Honghui Hospital Affiliated to Xi’an Jiaotong University. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and institutional requirements.
Author contributions
YS: Writing–original draft. GY: Writing–original draft. DL: Writing–original draft. YL: Formal Analysis, Writing–review and editing. CR: Formal Analysis, Writing–review and editing. YX: Data curation, Validation, Writing–review and editing. YY: Writing–review and editing. KZ: Writing–review and editing, Conceptualization. TM: Resources, Writing–review and editing. ZL: Funding acquisition, Resources, Writing–review and editing.
Funding
The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2023.1289976/full#supplementary-material
SUPPLEMENTARY TABLE S1 | Lists for the differentially expressed genes between OP and control samples in GSE56815 datasets.
SUPPLEMENTARY TABLE S2 | Results for the GSEA enrichment analysis of hub gene.
SUPPLEMENTARY TABLE S3 | Lists for the differentially expressed immune-related gene sets between OP and control samples in GSE56815 datasets.
SUPPLEMENTARY TABLE S4 | Results for the statistic analysis of spearman correlation analysis of differentially expressed immune-related gene sets (p value).
References
Aibar-Almazan A., Voltes-Martinez A., Castellote-Caballero Y., Afanador-Restrepo D. F., Carcelen-Fraile M. D. C., Lopez-Ruiz E. (2022). Current status of the diagnosis and management of osteoporosis. Int. J. Mol. Sci. 23 (16), 9465. doi:10.3390/ijms23169465
Cai P., Lu Y., Yin Z., Wang X., Zhou X., Li Z. (2021). Baicalein ameliorates osteoporosis via akt/foxo1 signaling. Aging (Albany NY) 13 (13), 17370–17379. doi:10.18632/aging.203227
Calvi L. M., Adams G. B., Weibrecht K. W., Weber J. M., Olson D. P., Knight M. C., et al. (2003). Osteoblastic cells regulate the haematopoietic stem cell niche. Nature 425 (6960), 841–846. doi:10.1038/nature02040
Che L., Wang Y., Sha D., Li G., Wei Z., Liu C., et al. (2023). A biomimetic and bioactive scaffold with intelligently pulsatile teriparatide delivery for local and systemic osteoporosis regeneration. Bioact. Mater. 19, 75–87. doi:10.1016/j.bioactmat.2022.03.023
Chen H., Shen G., Shang Q., Zhang P., Yu D., Yu X., et al. (2021). Plastrum testudinis extract suppresses osteoclast differentiation via the nf-κb signaling pathway and ameliorates senile osteoporosis. J. Ethnopharmacol. 276, 114195. doi:10.1016/j.jep.2021.114195
Chen L., Shi X., Weng S. J., Xie J., Tang J. H., Yan D. Y., et al. (2020). Vitamin K2 can rescue the dexamethasone-induced downregulation of osteoblast autophagy and mitophagy thereby restoring osteoblast function in vitro and in vivo. Front. Pharmacol. 11, 1209. doi:10.3389/fphar.2020.01209
Chen L., Wu X. Y., Jin Q., Chen G. Y., Ma X. (2023). The correlation between osteoporotic vertebrae fracture risk and bone mineral density measured by quantitative computed tomography and dual energy X-ray Absorptiometry: a systematic review and meta-analysis. Eur. Spine J. 32 (11), 3875–3884. doi:10.1007/s00586-023-07917-9
Cheng C. H., Chen L. R., Chen K. H. (2022). Osteoporosis due to hormone imbalance: an overview of the effects of estrogen deficiency and glucocorticoid overuse on bone turnover. Int. J. Mol. Sci. 23 (3), 1376. doi:10.3390/ijms23031376
Cheng L., Li H., Zhan H., Liu Y., Li X., Huang Y., et al. (2022). Alterations of M6a rna methylation regulators contribute to autophagy and immune infiltration in primary sjögren's syndrome. Front. Immunol. 13, 949206. doi:10.3389/fimmu.2022.949206
Cheng X., Yuan H., Cheng J., Weng X., Xu H., Gao J., et al. (2020). Chinese expert consensus on the diagnosis of osteoporosis by imaging and bone mineral density. Quantitative imaging Med. Surg. 10 (10), 2066–2077. doi:10.21037/qims-2020-16
Davis A. P., Grondin C. J., Johnson R. J., Sciaky D., Wiegers J., Wiegers T. C., et al. (2021). Comparative toxicogenomics database (ctd): update 2021. Nucleic Acids Res. 49 (D1), D1138–D1143. doi:10.1093/nar/gkaa891
Deng P., Yuan Q., Cheng Y., Li J., Liu Z., Liu Y., et al. (2021). Loss of Kdm4b exacerbates bone-fat imbalance and mesenchymal stromal cell exhaustion in skeletal aging. Cell Stem Cell 28 (6), 1057–1073.e7. doi:10.1016/j.stem.2021.01.010
Fan Y., Ge N., Wang X., Sun W., Mao R., Bu W., et al. (2014). Amplification and over-expression of Map3k3 gene in human breast cancer promotes formation and survival of breast cancer cells. J. Pathol. 232 (1), 75–86. doi:10.1002/path.4283
Fang E. F., Hou Y., Palikaras K., Adriaanse B. A., Kerr J. S., Yang B., et al. (2019). Mitophagy inhibits amyloid-beta and tau pathology and reverses cognitive deficits in models of alzheimer's disease. Nat. Neurosci. 22 (3), 401–412. doi:10.1038/s41593-018-0332-9
Filipowska J., Kondegowda N. G., Leon-Rivera N., Dhawan S., Vasavada R. C. (2022). Lgr4, a G Protein-Coupled receptor with a systemic role: from development to metabolic regulation. Front. Endocrinol. (Lausanne) 13, 867001. doi:10.3389/fendo.2022.867001
Filippiadis D. K., Charalampopoulos G., Mazioti A., Keramida K., Kelekis A. (2018). Bone and soft-tissue biopsies: what you need to know. Semin. Interv. Radiol. 35 (4), 215–220. doi:10.1055/s-0038-1669467
Fischer V., Haffner-Luntzer M. (2022). Interaction between bone and immune cells: implications for postmenopausal osteoporosis. Seminars cell & Dev. Biol. 123, 14–21. doi:10.1016/j.semcdb.2021.05.014
Foessl I., Dimai H. P., Obermayer-Pietsch B. (2023). Long-term and sequential treatment for osteoporosis. Nat. Rev. Endocrinol. 19 (9), 520–533. doi:10.1038/s41574-023-00866-9
Friedman J., Hastie T., Tibshirani R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33 (1), 1–22. doi:10.18637/jss.v033.i01
Golob A. L., Laya M. B. (2015). Osteoporosis: screening, prevention, and management. Med. Clin. North Am. 99 (3), 587–606. doi:10.1016/j.mcna.2015.01.010
Guo Y., Jia X., Cui Y., Song Y., Wang S., Geng Y., et al. (2021). Sirt3-Mediated mitophagy regulates ages-induced bmscs senescence and senile osteoporosis. Redox Biol. 41, 101915. doi:10.1016/j.redox.2021.101915
Han H., Xiao H., Wu Z., Liu L., Chen M., Gu H., et al. (2022). The mir-98-3p/jag1/notch1 Axis mediates the multigenerational inheritance of osteopenia caused by maternal dexamethasone exposure in female rat offspring. Exp. Mol. Med. 54 (3), 298–308. doi:10.1038/s12276-022-00743-x
Hanley J. A., McNeil B. J. (1982). The meaning and use of the area under a receiver operating characteristic (roc) curve. Radiology 143 (1), 29–36. doi:10.1148/radiology.143.1.7063747
He H., Cao S., Niu T., Zhou Y., Zhang L., Zeng Y., et al. (2016). Network-based meta-analyses of associations of multiple gene expression profiles with bone mineral density variations in women. PLoS One 11 (1), e0147475. doi:10.1371/journal.pone.0147475
He Y., Wang L., Liu W., Zhong J., Bai S., Wang Z., et al. (2015). Map3k3 expression in tumor cells and tumor-infiltrating lymphocytes is correlated with favorable patient survival in lung cancer. Sci. Rep. 5, 11471. doi:10.1038/srep11471
Hsu C. L., Lee W. C. (2010). Detecting differentially expressed genes in heterogeneous diseases using half student's T-test. Int. J. Epidemiol. 39 (6), 1597–1604. doi:10.1093/ije/dyq093
Hu Y., Tan L. J., Chen X. D., Liu Z., Min S. S., Zeng Q., et al. (2018). Identification of novel potentially pleiotropic variants associated with osteoporosis and obesity using the cfdr method. J. Clin. Endocrinol. Metab. 103 (1), 125–138. doi:10.1210/jc.2017-01531
Huang M., Ahmed A., Wang W., Wang X., Ma C., Jiang H., et al. (2022). Negative elongation factor (nelf) inhibits premature granulocytic development in zebrafish. Int. J. Mol. Sci. 23 (7), 3833. doi:10.3390/ijms23073833
Ito K., Murphy D. (2013). Application of Ggplot2 to pharmacometric graphics. CPT Pharmacometrics Syst. Pharmacol. 2 (10), e79. doi:10.1038/psp.2013.56
Jiang W. Y., Xing C., Wang H. W., Wang W., Chen S. Z., Ning L. F., et al. (2018). A lox/chop-10 crosstalk governs osteogenic and adipogenic cell fate by mscs. J. Cell Mol. Med. 22 (10), 5097–5108. doi:10.1111/jcmm.13798
Johnston C. B., Dagar M. (2020). Osteoporosis in older adults. Med. Clin. North Am. 104 (5), 873–884. doi:10.1016/j.mcna.2020.06.004
Kim H. L., Cho H. Y., Park I. Y., Choi J. M., Kim M., Jang H. J., et al. (2011). The positive association between peripheral blood cell counts and bone mineral density in postmenopausal women. Yonsei Med. J. 52 (5), 739–745. doi:10.3349/ymj.2011.52.5.739
Kim P. G., Niroula A., Shkolnik V., McConkey M., Lin A. E., Slabicki M., et al. (2021). Dnmt3a-Mutated clonal hematopoiesis promotes osteoporosis. J. Exp. Med. 218 (12), e20211872. doi:10.1084/jem.20211872
Krishnan A., Muthusami S. (2017). Hormonal alterations in pcos and its influence on bone metabolism. J. Endocrinol. 232 (2), R99–R113. doi:10.1530/JOE-16-0405
Kung A. W., Xiao S. M., Cherny S., Li G. H., Gao Y., Tso G., et al. (2010). Association of Jag1 with bone mineral density and osteoporotic fractures: a genome-wide association study and follow-up replication studies. Am. J. Hum. Genet. 86 (2), 229–239. doi:10.1016/j.ajhg.2009.12.014
Langdahl B. (2020). Treatment of postmenopausal osteoporosis with bone-forming and antiresorptive treatments: combined and sequential approaches. Bone 139, 115516. doi:10.1016/j.bone.2020.115516
Langfelder P., Horvath S. (2008). Wgcna: an R package for weighted correlation network analysis. BMC Bioinforma. 9, 559. doi:10.1186/1471-2105-9-559
Lee S. Y., An H. J., Kim J. M., Sung M. J., Kim D. K., Kim H. K., et al. (2021). Pink1 deficiency impairs osteoblast differentiation through aberrant mitochondrial homeostasis. Stem cell Res. Ther. 12 (1), 589. doi:10.1186/s13287-021-02656-4
Li D., Tian Y., Yin C., Huai Y., Zhao Y., Su P., et al. (2019). Silencing of lncrna Ak045490 promotes osteoblast differentiation and bone formation via beta-catenin/tcf1/runx2 signaling Axis. Int. J. Mol. Sci. 20 (24), 6229. doi:10.3390/ijms20246229
Li J., Yang S., Li X., Liu D., Wang Z., Guo J., et al. (2017). Role of endoplasmic reticulum stress in disuse osteoporosis. Bone 97, 2–14. doi:10.1016/j.bone.2016.12.009
Li L., Qi Y., Ma X., Xiong G., Wang L., Bao C. (2018). Trim22 knockdown suppresses chronic myeloid leukemia via inhibiting pi3k/akt/mtor signaling pathway. Cell Biol. Int. 42 (9), 1192–1199. doi:10.1002/cbin.10989
Liang G., Zhao J., Dou Y., Yang Y., Zhao D., Zhou Z., et al. (2022). Mechanism and experimental verification of luteolin for the treatment of osteoporosis based on network pharmacology. Front. Endocrinol. (Lausanne) 13, 866641. doi:10.3389/fendo.2022.866641
Liu F., Yuan L., Li L., Yang J., Liu J., Chen Y., et al. (2023). S-sulfhydration of Sirt3 combats bmsc senescence and ameliorates osteoporosis via stabilizing heterochromatic and mitochondrial homeostasis. Pharmacol. Res. 192, 106788. doi:10.1016/j.phrs.2023.106788
Liu T., Huang J., Xu D., Li Y. (2021). Identifying a possible new target for diagnosis and treatment of postmenopausal osteoporosis through bioinformatics and clinical sample analysis. Ann. Transl. Med. 9 (14), 1154. doi:10.21037/atm-21-3098
Liu Y., Li Y., Liu X., Wang C. S. (2019). Investigation of transcriptome mechanism associated with osteoporosis explored by microarray analysis. Exp. Ther. Med. 17 (5), 3459–3464. doi:10.3892/etm.2019.7349
Liu Y. C. G., Teng A. Y. (2023). Potential contribution of immature myeloid Cd11c(+)Dendritic cells-derived osteoclast precursor to inflammation-induced bone loss in the traf6-null chimeras in-vivo. J. Dent. Sci. 18 (3), 1372–1377. doi:10.1016/j.jds.2023.03.016
Liu Y. Z., Zhou Y., Zhang L., Li J., Tian Q., Zhang J. G., et al. (2015). Attenuated monocyte apoptosis, a new mechanism for osteoporosis suggested by a transcriptome-wide expression study of monocytes. PLoS One 10 (2), e0116792. doi:10.1371/journal.pone.0116792
Maduranga KWAH, Choi Y. H., Park S. R., Lee C. M., Kim G. Y. (2022). Bisphenol A inhibits osteogenic activity and causes bone resorption via the activation of retinoic acid-related orphan receptor α. J. Hazard Mater 438, 129458. doi:10.1016/j.jhazmat.2022.129458
Maity J., Barthels D., Sarkar J., Prateeksha P., Deb M., Rolph D., et al. (2022). Ferutinin induces osteoblast differentiation of dpscs via induction of Klf2 and autophagy/mitophagy. Cell Death Dis. 13 (5), 452. doi:10.1038/s41419-022-04903-9
Mattson M. P., Gleichmann M., Cheng A. (2008). Mitochondria in neuroplasticity and neurological disorders. Neuron 60 (5), 748–766. doi:10.1016/j.neuron.2008.10.010
Mizushima N., Levine B. (2020). Autophagy in human diseases. N. Engl. J. Med. 383 (16), 1564–1576. doi:10.1056/NEJMra2022774
Moayedi Y., Basch M. L., Pacheco N. L., Gao S. S., Wang R., Harrison W., et al. (2014). The candidate splicing factor sfswap regulates growth and patterning of inner ear sensory organs. PLoS Genet. 10 (1), e1004055. doi:10.1371/journal.pgen.1004055
Mohamad N. V., Ima-Nirwana S., Chin K. Y. (2021). The skeletal effects of gonadotropin-releasing hormone antagonists: a concise review. Endocr. Metab. Immune Disord. Drug Targets 21 (10), 1713–1720. doi:10.2174/1871530321666201216164410
Naik P. P., Birbrair A., Bhutia S. K. (2019). Mitophagy-driven metabolic switch reprograms stem cell fate. Cell Mol. Life Sci. 76 (1), 27–43. doi:10.1007/s00018-018-2922-9
Nogues X., Martinez-Laguna D. (2018). Update on osteoporosis treatment. Med. Clin. Barc. 150 (12), 479–486. doi:10.1016/j.medcli.2017.10.019
Park S. Y. (2018). Nomogram: an analogue tool to deliver digital knowledge. J. Thorac. Cardiovasc Surg. 155 (4), 1793. doi:10.1016/j.jtcvs.2017.12.107
Reyes J. M., Silva E., Chitwood J. L., Schoolcraft W. B., Krisher R. L., Ross P. J. (2017). Differing molecular response of young and advanced maternal age human oocytes to ivm. Hum. Reprod. 32 (11), 2199–2208. doi:10.1093/humrep/dex284
Salamanna F., Maglio M., Sartori M., Tschon M., Fini M. (2020). Platelet features and derivatives in osteoporosis: a rational and systematic review on the best evidence. Int. J. Mol. Sci. 21 (5), 1762. doi:10.3390/ijms21051762
Sarkar J., Das M., Howlader M. S. I., Prateeksha P., Barthels D., Das H. (2022). Epigallocatechin-3-Gallate inhibits osteoclastic differentiation by modulating mitophagy and mitochondrial functions. Cell Death Dis. 13 (10), 908. doi:10.1038/s41419-022-05343-1
Saxena Y., Routh S., Mukhopadhaya A. (2021). Immunoporosis: role of innate immune cells in osteoporosis. Front. Immunol. 12, 687037. doi:10.3389/fimmu.2021.687037
Scheibye-Knudsen M., Mitchell S. J., Fang E. F., Iyama T., Ward T., Wang J., et al. (2014). A high-fat diet and Nad(+) activate Sirt1 to rescue premature aging in cockayne syndrome. Cell Metab. 20 (5), 840–855. doi:10.1016/j.cmet.2014.10.005
Sebastian A., Contreras-Moreira B. (2013). The twilight zone of cis element alignments. Nucleic Acids Res. 41 (3), 1438–1449. doi:10.1093/nar/gks1301
Shannon P., Markiel A., Ozier O., Baliga N. S., Wang J. T., Ramage D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi:10.1101/gr.1239303
Shen W., Guo Y., Wang Y., Zhao K., Wang B., Yuille A. (2021). Deep differentiable random forests for age estimation. IEEE Trans. Pattern Anal. Mach. Intell. 43 (2), 404–419. doi:10.1109/TPAMI.2019.2937294
Sozen T., Ozisik L., Basaran N. C. (2017). An overview and management of osteoporosis. Eur. J. Rheumatol. 4 (1), 46–56. doi:10.5152/eurjrheum.2016.048
Sprangers S., Schoenmaker T., Cao Y., Everts V., de Vries T. J. (2016). Different blood-borne human osteoclast precursors respond in distinct ways to il-17a. J. Cell. physiology 231 (6), 1249–1260. doi:10.1002/jcp.25220
Srivastava R. K., Dar H. Y., Mishra P. K. (2018). Immunoporosis: immunology of osteoporosis-role of T cells. Front. Immunol. 9, 657. doi:10.3389/fimmu.2018.00657
Suarez-Farinas M., Lowes M. A., Zaba L. C., Krueger J. G. (2010). Evaluation of the psoriasis transcriptome across different studies by gene set enrichment analysis (gsea). PLoS One 5 (4), e10247. doi:10.1371/journal.pone.0010247
Tang C. Y., Wang H., Zhang Y., Wang Z., Zhu G., McVicar A., et al. (2022). GPR125 positively regulates osteoclastogenesis potentially through AKT-NF-κB and MAPK signaling pathways. Int. J. Biol. Sci. 18 (6), 2392–2405. doi:10.7150/ijbs.70620
Tao Z. S., Li T. L., Xu H. G., Yang M. (2021). Hydrogel contained valproic acid accelerates bone-defect repair via activating Notch signaling pathway in ovariectomized rats. J. Mater Sci. Mater Med. 33 (1), 4. doi:10.1007/s10856-021-06627-2
Unnanuntana A., Gladnick B. P., Donnelly E., Lane J. M. (2010). The assessment of fracture risk. J. Bone Jt. Surg. Am. 92 (3), 743–753. doi:10.2106/JBJS.I.00919
Wang H., Ou Y., Fan T., Zhao J., Kang M., Dong R., et al. (2021). Development and internal validation of a nomogram to predict mortality during the icu stay of thoracic fracture patients without neurological compromise: an analysis of the mimic-iii clinical database. Front. public health 9, 818439. doi:10.3389/fpubh.2021.818439
Wang C., Jiang X., Li Q., Zhang Y. Z., Tao J. F., Wu C. A. (2023). Identification of core pathogenic genes and pathways in elderly osteoporosis based on bioinformatics analysis. Zhonghua Yu Fang. Yi Xue Za Zhi 57 (7), 1040–1046. doi:10.3760/cma.j.cn112150-20230221-00140
Wang J., Zhang Y., Cao J., Wang Y., Anwar N., Zhang Z., et al. (2023). The role of autophagy in bone metabolism and clinical significance. Autophagy 19 (9), 2409–2427. doi:10.1080/15548627.2023.2186112
Wang S., Greenbaum J., Qiu C., Gong Y., Wang Z., Lin X., et al. (2023). Single-cell rna sequencing reveals in vivo osteoimmunology interactions between the immune and skeletal systems. Front. Endocrinol. (Lausanne) 14, 1107511. doi:10.3389/fendo.2023.1107511
Wang N., Wang Y., Zhang H., Guo Y., Chen C., Zhang W., et al. (2020). Association of bone mineral density with nine urinary personal care and consumer product Chemicals and metabolites: a national-representative, population-based study. Environ. Int. 142, 105865. doi:10.1016/j.envint.2020.105865
Wang S., Deng Z., Ma Y., Jin J., Qi F., Li S., et al. (2020). The role of autophagy and mitophagy in bone metabolic disorders. Int. J. Biol. Sci. 16 (14), 2675–2691. doi:10.7150/ijbs.46627
Wang Y., Liu L., Qu Z., Wang D., Huang W., Kong L., et al. (2022). Tanshinone ameliorates glucocorticoid-induced bone loss via activation of Akt1 signaling pathway. Front. Cell Dev. Biol. 10, 878433. doi:10.3389/fcell.2022.878433
Wang Z., Deng Z., Gan J., Zhou G., Shi T., Wang Z., et al. (2017). Tial(6)V(4) particles promote osteoclast formation via autophagy-mediated downregulation of interferon-beta in osteocytes. Acta Biomater. 48, 489–498. doi:10.1016/j.actbio.2016.11.020
Weng J., Yang Y., Song D., Huo R., Li H., Chen Y., et al. (2021). Somatic Map3k3 mutation defines a subclass of cerebral cavernous malformation. Am. J. Hum. Genet. 108 (5), 942–950. doi:10.1016/j.ajhg.2021.04.005
World Medical Association (2013). World medical association declaration of Helsinki: ethical principles for medical research involving human subjects. Jama 310 (20), 2191–2194. doi:10.1001/jama.2013.281053
Xia Q. D., Xun Y., Lu J. L., Lu Y. C., Yang Y. Y., Zhou P., et al. (2020). Network pharmacology and molecular docking analyses on lianhua qingwen capsule indicate Akt1 is a potential target to treat and prevent covid-19. Cell Prolif. 53 (12), e12949. doi:10.1111/cpr.12949
Xie L., Feng E., Li S., Chai H., Chen J., Li L., et al. (2023). Comparisons of gene expression between peripheral blood mononuclear cells and bone tissue in osteoporosis. Medicine 102 (20), e33829. doi:10.1097/md.0000000000033829
Xu C., Zhang Y., Wang Q., Xu Z., Jiang J., Gao Y., et al. (2016). Long non-coding rna Gas5 controls human embryonic stem cell self-renewal by maintaining nodal signalling. Nat. Commun. 7, 13287. doi:10.1038/ncomms13287
Xu W., Zhao D., Huang X., Zhang M., Yin M., Liu L., et al. (2022). The prognostic value and clinical significance of mitophagy-related genes in hepatocellular carcinoma. Front. Genet. 13, 917584. doi:10.3389/fgene.2022.917584
Xu Y., Shen J., Ran Z. (2020). Emerging views of mitophagy in immunity and autoimmune diseases. Autophagy 16 (1), 3–17. doi:10.1080/15548627.2019.1603547
Yan C., Shi Y., Yuan L., Lv D., Sun B., Wang J., et al. (2023). Mitochondrial quality control and its role in osteoporosis. Front. Endocrinol. (Lausanne) 14, 1077058. doi:10.3389/fendo.2023.1077058
Yang X., Jiang T., Wang Y., Guo L. (2019). The role and mechanism of Sirt1 in resveratrol-regulated osteoblast autophagy in osteoporosis rats. Sci. Rep. 9 (1), 18424. doi:10.1038/s41598-019-44766-3
Yin W., Shi L., Mao Y. (2019). Mir-194 regulates nasopharyngeal carcinoma progression by modulating Map3k3 expression. FEBS Open Bio 9 (1), 43–52. doi:10.1002/2211-5463.12545
Yu G., Wang L. G., Han Y., He Q. Y. (2012). Clusterprofiler: an R package for comparing biological themes among gene clusters. Omics 16 (5), 284–287. doi:10.1089/omi.2011.0118
Zanotti S., Canalis E. (2013). Notch signaling in skeletal health and disease. Eur. J. Endocrinol. 168 (6), R95–R103. doi:10.1530/EJE-13-0115
Zhang L., Dai L., Li D. (2021). Mitophagy in neurological disorders. J. neuroinflammation 18 (1), 297. doi:10.1186/s12974-021-02334-5
Zhang L., Chen X., Yin J. (2019). Prediction of potential mirna-disease associations through a novel unsupervised deep learning framework with variational autoencoder. Cells 8 (9), 1040. doi:10.3390/cells8091040
Zhang Y., Wang S. S., Tao L., Pang L. J., Zou H., Liang W. H., et al. (2019). Overexpression of MAP3K3 promotes tumour growth through activation of the NF-κB signalling pathway in ovarian carcinoma. Sci. Rep. 9 (1), 8401. doi:10.1038/s41598-019-44835-7
Zhao H., Wang D., Fu D., Xue L. (2015). Predicting the potential ankylosing spondylitis-related genes utilizing bioinformatics approaches. Rheumatol. Int. 35 (6), 973–979. doi:10.1007/s00296-014-3178-9
Zhao J., Lin F., Liang G., Han Y., Xu N., Pan J., et al. (2021). Exploration of the molecular mechanism of polygonati rhizoma in the treatment of osteoporosis based on network pharmacology and molecular docking. Front. Endocrinol. (Lausanne) 12, 815891. doi:10.3389/fendo.2021.815891
Zhao Y., Li L. (2023). The application of a topology optimization algorithm based on the kriging surrogate model in the mirror design and optimization of an aerial camera. Sensors (Basel) 23 (16), 7236. doi:10.3390/s23167236
Zong W. X., Rabinowitz J. D., White E. (2016). Mitochondria and cancer. Mol. Cell 61 (5), 667–676. doi:10.1016/j.molcel.2016.02.011
Keywords: bioinformatics, osteoporosis, mitophagy, differentially expressed genes (DEGs), biomarker, protein-protein interaction (PPI)
Citation: Su Y, Yu G, Li D, Lu Y, Ren C, Xu Y, Yang Y, Zhang K, Ma T and Li Z (2024) Identification of mitophagy-related biomarkers in human osteoporosis based on a machine learning model. Front. Physiol. 14:1289976. doi: 10.3389/fphys.2023.1289976
Received: 04 October 2023; Accepted: 21 December 2023;
Published: 08 January 2024.
Edited by:
Gang Wu, VU Amsterdam, NetherlandsReviewed by:
Zhiwen Luo, Fudan University, ChinaDengbao Yang, University of Texas Southwestern Medical Center, United States
Copyright © 2024 Su, Yu, Li, Lu, Ren, Xu, Yang, Zhang, Ma and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Teng Ma, gukemt@163.com; Zhong Li, sy1273341368@126.com
†These authors share first authorship