Identification of osteoarthritis-associated chondrocyte subpopulations and key gene-regulating drugs based on multi-omics analysis

Abstract

The mechanism by which chondrocytes respond to mechanical stress in joints significantly affects the balance and function of cartilage. This study aims to characterize osteoarthritis-associated chondrocyte subpopulations and key gene targets for regulatory drugs. To begin, single-cell and transcriptome datasets were obtained from the Gene Expression Omnibus (GEO) database. Cell communication and pseudo-temporal analysis, as well as High-dimensional Weighted Gene Co-expression Network Analysis (hdWGCNA), were conducted on the single-cell data to identify key chondrocyte subtypes and module genes. Subsequently, Consensus Cluster Plus analysis was utilized to identify distinct disease subgroups within the osteoarthritis (OA) training dataset based on the key module genes. Furthermore, differential gene expression analysis and GO/KEGG pathway enrichment analysis were performed on the identified subgroups. To screen for hub genes associated with OA, a combination of 10 machine learning algorithms and 113 algorithm compositions was integrated. Additionally, the immune and pathway scores of the training dataset samples were evaluated using the ESTIMATE, MCP-counter, and ssGSEA algorithms to establish the relationship between the hub genes and immune and pathways. Following this, a network depicting the interaction between the hub genes and transcription factors was constructed based on the Network Analyst database. Moreover, the hub genes were subjected to drug prediction and molecular docking using the RNAactDrug database and AutoDockTools. Finally, real-time fluorescence quantitative PCR (RT-qPCR) was employed to detect the expression of hub genes in the plasma samples collected from osteoarthritis patients and healthy adults. In the OA sample, there is a significant increase in the proportion of prehypertrophic chondrocytes (preHTC), particularly in subgroups 6, 7, and 9. We defined these subgroups as OA_PreHTC subgroups. The OA_PreHTC subgroup exhibits a higher communication intensity with proliferative-related pathways such as ANGPTL and TGF-β. Furthermore, two OA disease subgroups were identified in the training set samples. This led to the identification of 411 differentially expressed genes (DEGs) related to osteoarthritis, 2485 DEGs among subgroups, as well as 238 intersecting genes and 5 hub genes (MMP13, FAM26F, CHI3L1, TAC1, and CKS2). RT-qPCR results indicate significant differences in the expression levels of five hub genes and their related TFs in the clinical blood samples of OA patients compared to the healthy control group (NC). Moreover, these five hub genes are positively associated with inflammatory pathways such as TNF-α, JAK-STAT3, and inflammatory response, while being negatively associated with proliferation pathways like WNT and KRAS. Additionally, the five hub genes are positively associated with neutrophils, activated CD4 T cell, gamma delta T cell, and regulatory T cell, while being negatively associated with CD56dim natural killer cell and Type 17T helper cell. Molecular docking results reveal that CAY10603, Tenulin, T0901317, and Nonactin exhibit high binding activity to CHI3L1, suggesting their potential as therapeutic drugs for OA. The OA_PreHTC subgroups plays a crucial role in the occurrence and development of osteoarthritis (OA). Five hub genes may exert their effects on OA through interactions with PreHTC cells, other chondrocytes, and immune cells, playing a role in inhibiting cell proliferation and stimulating inflammation, thus having high diagnostic value for OA. Additionally, CAY10603, Tenulin, T0901317, and Nonactin have potential therapeutic effects for OA patients.

Keywords: Bioinformatics; Machine learning; Molecular docking; Osteoarthritis; Single cell analysis.

Fig. 1

Technology roadmap.

Fig. 2

Single-cell landscape analysis of osteoarthritis. (A) UMAP plots after cellular annotation of the osteoarthritis single-cell dataset (GSE152805); (B) Cellular proportions of the eight chondrocyte subpopulations in osteoarthritic and normal articular cartilage tissues; (C) Analysis of key marker genes and functional enrichment of the eight chondrocyte subpopulations.

Fig. 3

Classification of pre-HTC cells and pseudo-temporal analysis. (A) UMAP plot of cluster distribution difference between PreHTC in OA and Normal cartilage tissues; (B) Histogram of cluster distribution difference between PreHTC in OA and Normal cartilage tissues; (C) Pseudo-temporal differentiation trajectory map from Normal_PreHTC to OA_PreHTC; (D) Cell communication heatmap from Normal_PreHTC to OA_PreHTC.

Fig. 4

Cell communication analysis. (A) Cell communication results diagrams for OA_PreHTC and Normal_PreHTC as receptors; (B) Cell communication results diagrams for OA_PreHTC and Normal_PreHTC as ligands; (C) Cell communication network diagram for OA_PreHTC and Normal_PreHTC in the TGF-β signal; (D) Cell communication network diagram for OA_PreHTC and Normal_PreHTC in the ANGPTL signal; (E) Cell communication network diagram for OA_PreHTC and Normal_PreHTC in the SPP1 signal.

Fig. 5

Functional enrichment analysis of OA_PreHTC and Normal_PreHTC. A. Enrichment analysis of HYPOXIA, Angiogenesis, and Senescence signaling pathways in OA_PreHTC and Normal_PreHTC. B. SCENIC analysis in OA_PreHTC and Normal_PreHTC. C. Differential cell metabolism in OA_PreHTC and Normal_PreHTC. D. Enrichment analysis of Cellcall functions.

Fig. 6

The analysis of hdWGCNA in OA_PreHTC. (A) Soft threshold screen for hdWGCNA analysis; (B) Black and pink modules significantly enriched in 6, 7, and 9clusters; C. GO/KEGG enrichment analysis of black and pink module genes.

Fig. 7

Consensus Cluster Plus analysis. (A) Classification of OA diseased samples from the training set into cluster1 (n = 17) and cluster2 (n = 15) by consistent cluster analysis; (B) Differences in expression levels of black and pink module genes in clusters; (C) Volcano map of differentially expressed genes between clusters; (D) Distribution of differentially expressed genes between clusters; (E) GSVA pathway enrichment analysis between clusters.

Fig. 8

Screening hub genes through machine learning. (A) Validation efficacy of models evaluated based on different combinations of machine learning algorithms; (B) AUC curves for the training set; (C) AUC curves for the validation set GSE89408; (D) AUC curves for the validation set GSE82107; (E) AUC curves for the validation set GSE19060.

Fig. 9

The relationship between the hub genes and the immune and pathway. (A) Correlation of hub genes with pathways in the training set; (B) Correlation of hub genes with immune scores in the training set. (*in the figure represents P < 0.05, **represents P < 0.01, ***represents P < 0.001, ****represents P < 0.0001, and ns represents P > 0.05).

Fig. 10

Interaction network analysis of TF-hub genes. (A) The pseudo-temporal analysis of hub genes; (B) The regulatory network of hub genes and TFs.

Fig. 11

Interaction network analysis of hub genes-Drug. (A–E) Drug-hub gene regulatory networks were constructed for each of the five hub genes based on the RNAactDrug database. The selection parameters were FDR < 0.01, P-value < 0.001, and the data sources were CCLE, GDSC and Cellminer; (F–I) Molecular docking pattern map of CHI3L1 with drugs. The free energy of binding of CHI3L1 to CAY10603 was − 4.75 kcal/mol (F), the free energy of binding of CHI3L1 to Tenulin was − 6.55 kcal/mol (G), the free energy of binding of CHI3L1 to T0901317 was − 4.86 kcal/mol (H), the free energy of binding of CHI3L1 to The free energy of CHI3L1 binding to T0901317 was − 4.86 kcal/mol (H), and that of CHI3L1 binding to Nonactin was − 14.1 kcal/mol (I). In the figure, purple color represents CHI3L1, cyan color is the drug small molecule, green color is the amino acid residue, and red color is the hydrogen bond that binds them.

Fig. 12

The molecular expression levels of Hub genes and their associated TF were measured in clinical blood samples using RT-qPCR. (A–E) The mRNA expression levels of five hub genes (CHI3L1, CKS2, FAM26F, MMP13, and TAC1) in plasma of NC and OA groups; (F–I). The mRNA expression levels of hub genes-associated TFs (NFKB1, CREB1, CREB3, and RELA) in plasma of NC and OA groups; (J) The mRNA expression levels of HDAC6, a CAY10603 target, in plasma of NC and mRNA expression levels of HDAC6, a target of CAY10603, in plasma of OA group. (****P < 0.0001, vs. Control, N = 5, n = 3).