Machine learning algorithm-based biomarker exploration and validation of mitochondria-related diagnostic genes in osteoarthritis

Data collection

The study utilized data from the GEO, using GSE114007 as the experimental dataset to identify signature genes associated with OA, which included 18 normal and 20 OA cartilage samples. GSE57218 was used as the validation dataset to confirm the results of the first study, which included seven normal and 33 OA cartilage samples.

Identification of differentially expressed genes and mitochondria-related genes

Differentially expressed genes (DEGs) were analyzed with the R package “limma”. Screening criteria were —Log2FC— ¿= 1 and adjust P. Val < 0.05. A list of mitochondria-related genes was also obtained for subsequent studies. Specifically, the list of MRGs was obtained by searching for mitochondria in the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) using the keyword ‘mitochondria’, and the specific information is shown in Table S1.

Functional enrichment analysis of differentially expressed genes

To initially explore the functions and pathways played by DEGs in OA initiation and progression, enrichment analyses were performed. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “clusterProfiler” package. Adjust P.Val 0.05 indicates statistically significant differences. The results were visualized by using the “ggplot2” package, and the corresponding histograms and bubble plots were obtained.

Disease ontology analysis of differentially expressed genes

Disease ontology (DO) provides a unified description of genes from the perspective of diseases, which can be used to reveal the relationship and function between genes and diseases, thereby deepening the understanding of diseases and the development of therapeutic methods. In this study, disease ontology enrichment analysis of DEGs was performed using the “DOSE” package of the R, with the threshold set at adj.P.Val < 0.05. The analysis results were visualised using histograms and bubble plots.

LASSO regression analysis and SVM-RFE algorithm

To discover highly specific signature genes, this study combined two machine learning methods, LASSO and SVM-RFE, to improve the accuracy of signature gene selection. LASSO, implemented using the “glmnet” package, employs L1 regularization to effectively identify key features from high-dimensional linear data. SVM-RFE, implemented using the “e1071”, “kernlab”, and “caret” packages, recursively eliminates the least important features, making it suitable for handling nonlinear and small-sample data. By combining these methods, we leveraged their complementary strengths to enhance the robustness and precision of feature selection. Subsequently, a Venn diagram was used to identify the intersecting signature genes selected by both methods, ensuring the highest discriminatory signature genes were identified.

OA signature gene expression and diagnostic efficiency

To assess the accuracy and diagnostic value of the OA signature genes identified in Section ‘Differential expression analysis of signature genes and diagnostic efficacy for diseases’, we first analyzed the expression levels of the intersecting genes selected by LASSO and SVM-RFE algorithms using the “limma” and “ggpubr” packages and visualized the results using box plots. Following this, we employed the “pROC” package to generate receiver operating characteristic (ROC) curves and determine the area under the curve (AUC) values, thereby assessing the diagnostic efficacy of these signature genes for OA.

Immune infiltration analysis

It is worth noting immune cell-mediated inflammatory responses significantly contribute to the development of OA. Single-sample gene set enrichment analysis (ssGSEA) was conducted to preliminarily investigate the role of immune cells in the disease process of OA. In this study, we used the ‘GSVA’ package to conduct ssGSEA analysis to clarify the immune landscape in OA cartilage tissues and the relationship between signature gene expression and the level of immune cell infiltration. Twenty-eight types of immune cells and their characteristic markers were identified from the results of previous studies (Subramanian et al., 2005). Beyond the initial analysis with the GSE114007 dataset, we also conducted immune cell infiltration and GRPEL1 correlation analysis on the validation dataset GSE57218.

Cartilage primary cell culture and cell models of OA

The isolation of mouse chondrocytes was performed according to previously published methods (Jin et al., 2021). Specifically, 10 four-week-old healthy and normal C57BL/6J mice were obtained from Zhengzhou University and housed in the SPF animal house. The mice were euthanized at the time of the experiment using an overdose of sodium pentobarbital. Samples of articular cartilage from the knee and femoral heads of mice were collected, cut into pieces, and washed with sterile PBS. This was followed by digestion using trypsin for 30 min, removal of trypsin, and washing with PBS. The samples were then digested using type II collagenase for 5 h, filtered, and the liquid retained. Finally, culture medium was added for incubation. Chondrocytes were cultured in DMEM/F12 medium (10091; Procell, Wuhan, China) containing 10% FBS (10091; Gibco, Billings, MT, USA) and maintained at 37 °C with 5% CO₂. To simulate osteoarthritis in vitro, the second-generation chondrocytes were treated with 5 μg/mL LPS (L2880-10MG; Sigma-Aldrich, St. Louis, MO, USA) based on literature references (Chen et al., 2020; Li & Chen, 2021). The cells were treated with LPS for 24 h as the experimental group, while the untreated cells were used as the control group. The success of the model was confirmed by RT-qPCR, and subsequent experiments were conducted using third-generation cells. The guidelines for the care and use of laboratory animals are conducted in all procedures involving animals. Measures to minimize animal suffering included the use of anesthesia during chondrocyte isolation and providing appropriate post-operative care. And the animal experiment protocol was approved by the Zhengzhou Weisha Biotechnology Co. Laboratory Animal Ethics Committee (Approval number: V3A02022109012).

Quantitative Real-time PCR

Total RNA was extracted from cells using the Total RNA Kit I (R6834-02; Omega, Norcross, GA, USA), and its quality and purity were assessed via NanoDrop (Thermo Fisher Scientific, Waltham, MA, USA). The extracted RNA was then reverse-transcribed into cDNA using the NovoScript Plus cDNA SuperMix kit (E047-01B; Novoprotein, Beijing, China). The expression of the signature gene was quantified using the Novostart SYBR qPCR SuperMix (E096-01A; Novoprotein), and the 2^−ΔΔCT method was employed for result calculation. The primer sequence can be found in Table S2.

Cell Transfection

One day before transfection, OA cells were seeded in 6-well plates to reach 70–80% confluency at the time of transfection. Transfection was performed using Lipofectamine 3000 (L3000015; Invitrogen) according to the manufacturer’s instructions. The experiment was divided into two groups: the control group (NC group, pCMV-EGFP- 3 ×FLAG-Neo vector) and the overexpression group (OE group, pCMV-EGFP-GRPEL1- 3 ×FLAG-Neo vector). For each well, 2.5 µg of the respective plasmid and 5 µl of Lipofectamine 3000 were mixed and added. After 6 h of transfection, the medium was replaced with complete medium containing 10% FBS, and the cells were cultured for an additional 48 h. Transfection efficiency was verified by RT-PCR to measure GRPEL1 expression levels.

Statistical methods

The statistical analyses conducted in this study utilized the R language (version 4.1.1; R Core Team, 2021). A comparison between the two groups was made using either the student’s t-test or the Wilcoxon test. Spearman’s correlation was employed to assess the relationship between the two variables. A significance level of p < 0.05 was deemed indicative of statistical significance.

Results of differentially expressed genes

To explore the genes that play potentially key roles in the pathological process of OA, the data microarray set GSE114007 was first subjected to differential expression analysis, and the results revealed a vast number of genes with abnormal expression, with a total of 1,502 DEGs in the OA cartilage samples. Among them, there were 790 genes with increased expression and 712 genes with decreased expression in OA, as shown in Fig. 1 and Table S1.

Functional enrichment analysis of differentially expressed genes

We conducted a GO analysis to annotate their functions based on the above-screened DEGs. In this analysis, we included all differentially expressed genes, both upregulated and downregulated, to preliminarily explore their biological functions in the pathological process of OA. Among them, biological process was mainly enriched in extracellular matrix organisation, extracellular structure organisation, external encapsulating strcture organisation, ossification, etc. Cellular component is mainly enriched in collagen-containing extracellular matrix, basement membrane, collagen trimer, endoplasmic reticulum lumen, transcription regulator complex, etc. Molecular function is mainly enriched in extracellular matrix structural constituent, conferring tensile strength glycosaminoglycan binding, heparin binding, DNA-binding transcription activator activity, RNA-polymerase II-specific, sulfur compound binding, and so on (Figs. 2A and 2B). These enriched processes are closely related to the pathological progression of OA. The degradation and remodeling of the extracellular matrix are hallmark features of OA, while abnormal ossification accelerates cartilage degeneration (Rahmati et al., 2017; Fujii et al., 2022). Additionally, the degradation of collagen and glycosaminoglycans significantly disrupts the articular cartilage’s structural integrity and mechanical properties. This degradation process leads to the loss of cartilage matrix and functional impairment, resulting in joint instability, uneven load distribution, and cartilage erosion, ultimately driving the pathological progression of OA (Charlier et al., 2019; Chen et al., 2023). These findings provide crucial biological insights into the pathological mechanisms of OA.

After that, to explore the pathways through which DEGs were involved in the pathological process of OA, KEGG analysis was also conducted in this study, and the results are shown in Figs. 2C and 2D. Specifically, they were mainly enriched in signaling pathways such as the FoxO signaling pathway, ECM-receptor interaction, PI3K-Akt signaling pathway, HIF-1 signaling pathway, and TNF signaling pathway.

Disease ontology analysis of differentially expressed genes

To understand whether DEGs may function in the pathological process of OA, DO analysis was performed on the DEGs screened above. The results suggested DEGs were mainly enriched in diseases such as osteoarthritis, tumors of the musculoskeletal system, and malignant tumors of other systems (Fig. 3). This also suggests the presence of potential signature genes of OA in these DEGs.

LASSO regression analysis and SVM-RFE algorithm to select OA signature genes

The aforementioned findings suggest the presence of certain individuals within the differentially expressed genes (DEGs) who possess crucial roles in the initiation and advancement of OA. Consequently, this study aims to conduct comprehensive analyses to identify the signature genes associated with OA. Firstly, the results of differential expression analysis revealed 1502 DEGs in the cartilage tissues of OA patients. Secondly, to further narrow down the range of feature genes, the intersection of DEGs and MRGs was taken, and 45 potential signature genes for the disease were finally obtained (Fig. 4A). And the specific list of signature genes is shown in Table S3. Additionally, enrichment analysis of MRGs within the DEGs revealed significant involvement in mitochondrial transmembrane transport, regulation of mitochondrial organization, and oxidoreductase activity. These findings suggest that mitochondrial dysfunction plays a pivotal role in the progression of OA, providing potential targets for therapeutic intervention (Fig. S1). Subsequently, 10 and 2 potential signature genes were determined respectively by using LASSO analysis and SVM-RFE algorithms (Figs. 4B, 4C). Finally, after cross-referencing the above results, two genes (MTFP1, GRPEL1) were identified. Therefore, these 2 genes were initially identified as signature genes of OA (Fig. 4D).

Differential expression analysis of signature genes and diagnostic efficacy for diseases

To specify the expression levels of OA signature genes in cartilage tissues, differential expression analysis was conducted. Firstly, the expression of MTFP1 and GRPEL1 were extracted from the experimental group dataset GSE114007 for differential analysis. The results revealed the expression of MTFP1 and GRPEL1 were abnormally downregulated in the OA group (Figs. 5A and 5B). However, it is not known whether they are specific and whether they have diagnostic value for OA. Therefore, ROC curves were plotted to clarify this issue. However, it is not known whether they are specific and whether they have diagnostic value for OA. Therefore, ROC curves were plotted to clarify this issue. The results showed AUC of GRPEL1 was 0.994, and that of MTFP1 was 1 (Figs. 5C, 5D). This implies that GRPEL1 and MTFP1 have good diagnostic values for OA. In conclusion, the above results suggest GRPEL1 and MTFP1 can be used as signature genes of OA and have better diagnostic efficacy.

Validation of signature gene expression and diagnostic efficacy

To clarify the accuracy of the above findings, the expression of MTFP1 and GRPEL1 in the test group dataset GSE57218 was extracted for differential analysis. The results showed only GRPEL1 expression was abnormally down-regulated in the OA group (see Figs. 6A, 6B). Similarly, the diagnostic efficacy of GRPEL1 for OA was further verified by ROC curves, and the AUC value was determined to be 0.797 (Figs. 6C, 6D). This indicated that GRPEL1 also had some diagnostic value for OA in the test group.

Further, to further verify the reliability of the results, this study used mouse chondrocytes to construct an in vitro OA model to verify the expression of the signature genes. First, the changes of extracellular matrix and degradative enzymes in chondrocytes after 5 ug/ml LPS treatment were detected by RT-qPCR, and the results revealed the expression of Ifngr2, Adamts5 and Mmp13 increased, and the expression of aggrecan and collagen II decreased, suggesting the model was successfully constructed, as shown in Fig. 7A. Subsequently, RT-qPCR was performed to detect the expression of Grpel1, and the results indicated its expression was reduced in the experimental group, as shown in Fig. 7B.

The above results suggest that GRPEL1 has a reliable diagnostic efficacy for OA. Specifically, GRPEL1 expression level was found to be a useful diagnostic marker based on our analysis.

Impact of Grpel1 overexpression on chondrocyte function

To further investigate the potential role of Grpel1, primary OA cells were transfected to overexpress Grpel1, and the expression changes of representative genes, including extracellular matrix components, matrix-degrading enzymes, and inflammatory molecules, were detected using RT-qPCR. These genes included Mmp13, Adamts5, Col2a1, Ifngr2, and Aggrecan. The results are shown in Fig. 8. Compared to the control group, Grpel1 overexpression significantly upregulated the expression of Col2a1 and Aggrecan, while significantly downregulating the expression of Mmp13, Adamts5, and Ifngr2. These findings suggest that Grpel1 may play a protective role in the pathogenesis of OA by promoting cartilage matrix synthesis, inhibiting matrix degradation, and reducing inflammatory responses.

Relationship between signature gene expression and level of immune cell infiltration

Immune cells are important in the prevention and treatment of OA through the secretion of cytokines and chemokines that promote inflammatory responses and articular cartilage destruction (Chen et al., 2021a). To investigate the degree of immune cell infiltration in OA, ssGSEA analysis was conducted on the experimental group data set. First, a difference analysis was undertaken on the level of immune cell infiltration in the OA and control groups, as shown in Fig. 9A. The results revealed the infiltration levels of activated B cells, eosinophils, type 17 T helper cells, memory B cells, central memory CD8 T cells, and effector memory CD8 T cells were decreased in OA samples compared to normal samples. In contrast, infiltration levels of activated dendritic cells, gamma-delta T cells, immature dendritic cells, MDSCs, macrophages, natural killer T cells, regulatory T cells, T follicular helper cells, and type 2 T helper cells were increased in OA samples. Summarily, multiple immune cells undergo infiltration abundance alterations in the OA immune landscape. To initially clarify whether the involvement of GRPEL1 in the pathological process of OA is related to immune cells, further correlation analyses were undertaken. The findings of this study demonstrate a positive correlation between GRPEL1 and the infiltration level of certain immune cells, namely MDC and macrophages, while a negative correlation was observed with the infiltration level of NK cells, as depicted in Fig. 9B.

Additionally, immune cell infiltration and correlation analyses were conducted based on the validation dataset GSE57218, as shown in Fig. S2. Specifically, the infiltration abundance of MDSCs, macrophages, and T follicular helper cells was higher in OA tissues, consistent with the results from the experimental dataset. However, the correlation analysis between GRPEL1 and immune cell infiltration did not show consistency with the results from the experimental dataset.