Transcriptomic analyses and machine-learning methods reveal dysregulated key genes and potential pathogenesis in human osteoarthritic cartilage

Abstract

Aims: This study aimed to explore the biological and clinical importance of dysregulated key genes in osteoarthritis (OA) patients at the cartilage level to find potential biomarkers and targets for diagnosing and treating OA.

Methods: Six sets of gene expression profiles were obtained from the Gene Expression Omnibus database. Differential expression analysis, weighted gene coexpression network analysis (WGCNA), and multiple machine-learning algorithms were used to screen crucial genes in osteoarthritic cartilage, and genome enrichment and functional annotation analyses were used to decipher the related categories of gene function. Single-sample gene set enrichment analysis was performed to analyze immune cell infiltration. Correlation analysis was used to explore the relationship among the hub genes and immune cells, as well as markers related to articular cartilage degradation and bone mineralization.

Results: A total of 46 genes were obtained from the intersection of significantly upregulated genes in osteoarthritic cartilage and the key module genes screened by WGCNA. Functional annotation analysis revealed that these genes were closely related to pathological responses associated with OA, such as inflammation and immunity. Four key dysregulated genes (cartilage acidic protein 1 (CRTAC1), iodothyronine deiodinase 2 (DIO2), angiopoietin-related protein 2 (ANGPTL2), and MAGE family member D1 (MAGED1)) were identified after using machine-learning algorithms. These genes had high diagnostic value in both the training cohort and external validation cohort (receiver operating characteristic > 0.8). The upregulated expression of these hub genes in osteoarthritic cartilage signified higher levels of immune infiltration as well as the expression of metalloproteinases and mineralization markers, suggesting harmful biological alterations and indicating that these hub genes play an important role in the pathogenesis of OA. A competing endogenous RNA network was constructed to reveal the underlying post-transcriptional regulatory mechanisms.

Conclusion: The current study explores and validates a dysregulated key gene set in osteoarthritic cartilage that is capable of accurately diagnosing OA and characterizing the biological alterations in osteoarthritic cartilage; this may become a promising indicator in clinical decision-making. This study indicates that dysregulated key genes play an important role in the development and progression of OA, and may be potential therapeutic targets.

Fig. 1

Identification of overlapping characteristic genes. a) Volcano plot of differentially expressed genes (DEGs) between osteoarthritis (OA) samples and control samples (logFC > 1 and adjusted p-value < 0.05). b) Soft-threshold power for weighted gene coexpression network analysis (WGCNA). The red line in the left panel indicates scale-free topological fit index = 9. c) Clustering dendrograms of all expressed genes with dissimilarity based on the topological overlap along with the different assigned modules. d) Heatmap of the correlations between modules and clinical traits. Red represents positive correlations, and blue represents negative correlations. e) A scatterplot of gene significance versus module membership in the green-yellow module. f) Venn diagram of the intersection of the overlapping characteristic genes obtained from upregulated DEGs and key module genes. g) Heatmap of the overlapping characteristic genes.

Fig. 2

Functional annotation of the overlapping characteristic genes. a) Gene Ontology functional analysis of the overlapping characteristic genes. b) Kyoto Encyclopedia of Genes and Genomes pathway analysis of the overlapping characteristic genes.

Fig. 3

Identification of key characteristic genes. a) Characteristic genes that were closely related to the study group were screened with the Boruta algorithm (green and yellow). b) Validation of the optimal gene expression signature by support vector machine–recursive feature elimination (SVM–RFE) algorithm selection. The colours in the figure represent the corresponding model accuracy for different gene numbers. c) The optimal lambda value was determined when the misclassification reached the minimum through the lasso regression algorithm. d) Venn diagram of the key characteristic genes intersected by the SVM–RFE and LASSO algorithms.

Fig. 4

Diagnostic value of the key characteristic genes. a) Receiver operating characteristic (ROC) analysis of the independent diagnostic efficacy of the four key genes in the training cohort (GSE169077, GSE57218). b) ROC analysis of the overall diagnostic efficacy of the four key genes in the training cohort. c) ROC analysis of the independent diagnostic efficacy of the four key genes in the external validation cohort (GSE89408, GSE143514). d) ROC analysis of the overall diagnostic efficacy of the four key genes in the external validation cohort. e) Decision curve analysis to evaluate the potential clinical value in the training cohort. f) Decision curve analysis to evaluate the potential clinical value in the external validation cohort. g) Differential expression levels of the four key genes between the osteoarthritis (OA) and normal groups in the external validation cohort (Mann–Whitney U test).

Fig. 5

Potential effect of key genes in osteoarthritis (OA). a) Coexpression network of the characterized genes. b) Gene Ontology functional (GO) analysis of the coexpressed genes. c) Kyoto Encyclopedia of Genes and Genomes pathway analysis (KEGG) of the coexpressed genes. d) Boxplots of the risk score between OA and healthy controls in the training cohort (Mann-Whitney U test).

Fig. 6

e) Boxplots of the risk score between osteoarthritis (OA) and healthy controls in the external validation cohort (Mann–Whitney U test). f) to i) Gene set enrichment analysis conducted between high- and low-risk groups. Specific f) Gene Ontology (GO) biological processes and g) Kyoto encyclopedia of Genes and Genomes (KEGG) pathways in the high-risk group of the training cohort (GSE169077, GSE57218). Specific h) GO biological processes and i) KEGG pathways in the high-risk group of the external validation cohort (GSE89408, GSE143514).

Fig. 7

Correlations of key genes with pathological changes in osteoarthritis (OA). a) Correlation analysis among the four key genes; darker red indicates stronger positive correlations. b) The correlations of four key genes with matrix metalloproteinases (MMPs), a disintegrin and metalloproteinase with thrombospondin motifs 5 (ADAMTS-5), ACAN, alkaline phosphatase (ALPL), and collagen type I alpha 1 (COL1A1); darker red indicates stronger positive correlations and darker blue indicates stronger negative correlations. c) The specific value of the correlation and statistical significance of the four key genes with MMPs, ADAMTS-5, ACAN, ALPL, and COL1A1 (Spearman correlation analysis; *p < 0.05, **p < 0.01, and ***p < 0.001).

Fig. 8

The landscape of immune infiltration in association with osteoarthritis (OA). a) Violin plot exhibiting the estimated immune score between OA and healthy controls (Mann–Whitney U test). b) Boxplots of the different immune cell infiltration profiles of OA and healthy controls (Mann–Whitney U test). c) The relationship between the four key genes and immune cell abundance (Spearman correlation analysis; *p < 0.05, **p < 0.01, and ***p < 0.001).