Identification of inflammation‑related biomarkers for osteoarthritis diagnosis and stratification through bulk and single‑cell RNA‑sequencing integration

Introduction

Osteoarthritis (OA) is a common chronic degenerative joint disease that primarily affects adults, with incidence increasing with age and affecting >500 million individuals worldwide (1,2). The pathogenesis of OA is yet to be elucidated and consequently, no curative therapy is available. However, although a number of treatment options are available, postoperative rehabilitation following total joint arthroplasty poses persistent challenges, including muscle weakness and atrophy caused by prolonged immobility, which impair joint stability and function, as well as joint stiffness and reduced range of motion resulting from scar tissue formation or disuse (3). Among pharmacological options, oral and topical non-steroidal anti-inflammatory drugs (NSAID), including cyclooxygenase-2 inhibitors, are effective for alleviating pain but carry increased risks. NSAIDs may cause gastrointestinal ulcers, bleeding or perforation, particularly with long-term use or in patients with a history of peptic ulcer disease. NSAIDs may increase the risk of cardiovascular complications, including hypertension, heart failure exacerbation, myocardial infarction and stroke (4). Biomaterial-based targeted delivery systems release drugs in response to specific internal or external stimuli associated with the disease microenvironment, but an absence of standardized safety and biocompatibility guidelines limits their use in the general population. For example, spherical lipid bilayer vesicles are capable of encapsulating both hydrophilic and hydrophobic drugs, thereby improving their stability and facilitating targeted delivery (5). Targeted strategies have focused on regulating pathways such as the PI3K/AKT/mTOR, EGFR, toll-like receptors (TLRs) and proinflammatory IL pathways in cartilage (6-8). Notably, intra-articular therapy holds potential for OA treatment as numerous cDNA-based interventions exhibit therapeutic effects in mouse models (9). Therefore, screening for effective biomarkers is important for guiding therapeutic development. Increasing evidence indicates the inflammatory nature of OA (10,11). Inflammation serves a key role in OA pathogenesis. Inflammatory mediators, such as IL-1b, IL-6 and TNF-a, derived from synovial tissues, cartilage and subchondral bone contribute to OA development and increase with age (12), which is consistent with the increased prevalence of OA in adults. Synovitis is commonly observed in OA-affected joints and is associated with radiographic progression and pain severity (13,14). In zymosan-induced arthritis, low-density lipoprotein-1 (LOX-1)-knockout mice exhibit decreases in inflammation, synovial hyperplasia, cartilage degeneration and MMP-3 expression levels (15). In a joint destabilization model of OA, LOX-1-knockout mice have attenuated disease with decreased levels of oxidized LDL (ox-LDL), runt-related transcription factor 2 (Runx2) and type X collagen, implicating LOX-1/ox-LDL-driven endochondral ossification (16). Cartilage inflammation is also associated with chloride-channel dysfunction (17). One study has revealed the involvement of pathophysiological inflammatory biomarkers in OA disease, including MMPs and TIMPs (18). Therefore, the present study suggested that inflammation-related biomarkers may represent possible avenues for OA diagnosis and treatment.

Bioinformatics is an important tool for integrating and interrogating biological data, which enables researchers to investigate novel research directions (19). Bulk RNA sequencing (RNA-seq) profiles average gene expression levels across mixed cell populations, enabling differential expression analysis but masking cell-type-specific signals (20). However, single-cell (sc)RNA-seq enables the characterization of functional cellular states of different cell types, cellular plasticity upon stimulation and cell differentiation or reprogramming trajectories (21). A number of previous studies use bulk RNA-seq (22-24) or scRNA-seq alone (25-27) to investigate the mechanisms underlying OA development and progression. However, integrating bulk and scRNA-seq provides a more comprehensive overview. Therefore, the present study combined these approaches in order to elucidate the molecular and cellular mechanisms driving OA progression. Furthermore, four machine learning (ML) models, namely, generalized linear models (GLM), random forests (RF), support vector machines (SVM) and extreme gradient boosting (XGB), were used to construct predictive models and possibly elucidate complex biological mechanisms. Therefore, the present study aimed to integrate bioinformatics and ML algorithms in order to highlight possible OA biomarkers and inflammation-related genes, followed by a scRNA-seq-based annotation of these biomarkers.

Materials and methods

Data collection

All human OA transcriptomic datasets, including bulk and scRNA-seq data, were obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). Transcriptome samples from the OA and control groups were screened using the search terms ‘osteoarthritis’ and ‘human’, respectively. The primary analysis incorporated six datasets, GSE12021(28), GSE1919(29), GSE36700(30), GSE51588(31), GSE55235(32) and GSE55457(32), comprising 99 samples. Although GSE12021 and GSE1919 were generated primarily for rheumatoid arthritis (RA) synovial studies, they also contain data on non-RA control samples, including OA and normal tissues. In the present study, OA and normal samples were selected and no RA samples were included in the present analyses. The validation group consisted of three datasets, GSE82107(33), GSE117999(34) and GSE169077(35), containing a total of 52 samples. The scRNA-seq dataset used in the present study was GSE104782(36), which was sourced from the GPL20301 platform (Illumina, Inc.). The patient clinical information is summarized in Table SI.

Identification of inflammation-related differentially expressed genes (DEGs)

To integrate OA and control samples from the six datasets, the ‘ComBat’ function in the ‘sva’ R package (37) (version 3.52.0; Posit Software, PBC) was used to eliminate batch differences in gene expression across datasets. Briefly, ‘ComBat’ was applied with a batch set to the dataset source, included no additional covariates in the model in order to preserve biological variance and used the parametric empirical Bayes framework (38). The processed dataset contained 60 OA samples and 39 control samples. DEGs were identified using an absolute log2-fold change ≥1 and a false discovery rate <0.05. To characterize inflammation-related biology among these DEGs, the curated ‘inflammatory response’ gene set was retrieved from the ‘c2.cp.kegg.symbols.gmt’ collection in the Molecular Signatures Database (version 2023.2; Broad Institute, Inc.) and these genes were intersected with DEGs to obtain inflammation-related DEGs (IRDGs). The chromosomal locations of the intersecting genes and their protein-protein interaction (PPI) networks were mapped and analyzed.

Cluster analysis and immune characterization

On the basis of the expression profiles of the IRDGs, a consistency clustering analysis (39) was carried out using the ‘ConsensusClusterPlus’ R/Bioconductor package (version 1.68.0; Matt Wilkerson and Peter Waltman) (40,41). Cluster separation was visualized using heatmaps and assessed through principal component analysis (PCA) (42). Gene set variation analysis (GSVA) (43) was carried out to assess the pathway variability among different clusters. The cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) R script (44) was used to deconvolute the gene expression matrix into the relative proportions of immune cell types. The estimation of stromal and immune cells in malignant tumor tissues using expression data R package (version 1.0.13; R Core Team) (45) was used to calculate immune infiltration scores (immune score), overall stromal content score (stromal scores) and composite scores (estimate scores) on the basis of gene expression data in OA samples (46). Pearson's correlation coefficients were used to assess the associations between these scores and different types of immune cell populations (47). Differences in immune cell distribution, immune function and immune checkpoint expression among the different clusters were characterized based on immune infiltration results.

ML model construction and validation

Initial screening of IRDGs was carried out using least absolute shrinkage and selection operator (LASSO) regression with the ‘glmnet()’ function in the ‘glmnet’ R package (48) (version 4.1.8; CRAN; https://cran.r-project.org/package=glmnet). The penalty parameter was set to α=1 (L1 regularization). The optimal λ was determined by tenfold cross-validation and the final model used ‘λ_min’, which was the value yielding the lowest cross-validation error. The selected genes were used as inputs for four machine learning models, namely GLM, RF, SVM and XGB. After data partitioning, a prediction function was constructed to analyze disease-specific signature genes. Predictors were standardized using the ‘glmnet’ default settings (standardize=TRUE). Subsequently, the following models and parameters were used: SVM [‘e1071’ version 1.7-16, C=1.0, γ=‘scale’ and 5-fold cross validation (CV)]; GLM/LASSO (‘glmnet’ version 4.1-8, family=‘binomial’ and α=1); RF (version 4.7-1.1, repeated CV with 5 folds and savePredictions=TRUE); and XGB (‘xgboost’ version 1.7.8.1, booster=‘gbtree’, max_depth=6 and learning_rate=0.1). The predictive accuracies of the four algorithms were compared by evaluating the residual box plots, reverse cumulative distribution of residuals and receiver operating characteristic (ROC) curves (49). The four most important genes were identified on the basis of their gene importance scores. Calibration curves and nomograms (50) were constructed for the best-performing models. To validate the prediction accuracy, the final model was validated using three independent datasets, GSE82107, GSE117999 and GSE169077, in which each inflammation-related gene was used as an individual predictor within the finalized model.

scRNA-seq data analysis

Single-cell gene attainment matrices were converted into Seurat object data structures using the ‘Seurat’ R package (version 4.4.0; CRAN). After integration and batch correction, gene expression was normalized for each cell using the Seurat ‘LogNormalize’ method (scaling factor=10,000). The Seurat ‘ScaleData’ was then applied to center and scale gene expression across all genes to support downstream analyses. Samples containing <3 cells and <50 genes were removed. Samples that met these thresholds were retained unless they failed additional quality control (QC) criteria. Strict quality control was applied to exclude low-quality cells, which were defined as cells with <200 detected genes or mitochondrial gene expression accounting for 10% or more of total transcripts (51). The data were batch corrected and normalized, and the top 3,000 genes with the highest intercellular coefficients of variation were extracted for further analysis. Following PCA and dimensionality reduction, t-distributed stochastic neighbor embedding (t-SNE) (52) was used to calculate the difference in neighbor distance for cell clustering. Marker genes for each subpopulation were identified and disease-specific signature gene expression was visualized. The ‘SingleR’ R package (version 2.6.0) (53) was used to annotate cell types. Single-cell differentiation trajectory analysis was carried out using the ‘monocle’ R package (54) to generate trajectory maps of dendritic cells, temporal changes, cell-type annotation and clustering (55).

Cell culture and transfection

ATDC5 cells (The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences) were cultured in a 1:1 mixture of DMEM and Ham's F-12 (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS (FuHeng Biology Co., Ltd.), 100 U/ml penicillin and 100 µg/ml streptomycin, and incubated at 37˚C in a humidified 5% CO2 environment. Lentiviral short hairpin RNAs (shRNAs) targeting low-density lipoprotein receptor (LDLR), adrenomedullin (ADM), MYC and NF-κB inhibitor-α (NFKBIA) as well as a non-targeting shRNA control were purchased from Shanghai GenePharma Co., Ltd. In addition, shRNA oligonucleotides were cloned into the pLKO.1-puro lentiviral backbone (Sigma-Aldrich; Merck KGaA). The following shRNA sequences (in the 5'-3' direction) were used: shRNA control, TACTGTTAGGATCAGGAGG; LDLR shRNA#1, CGGAGGTGACCAACAATAGAA; LDLR shRNA#2, GATGTCATAAACGAAGCCATT; ADM shRNA#1, CCTCATTACTACTTGAACTTT; ADM shRNA#2, TCCAGACTCTTTAGGATATAG; MYC shRNA#1, CGACGAGGAAGAGAATTTCTA; MYC shRNA#2, GCTTCGAAACTCTGGTGCATA; NFKBIA shRNA#1, GAGTCAGAATTCACAGAGGAT; and NFKBIA shRNA#2, GGACGAGGAGTACGAGCAAAT. Lentiviral particles were produced in 293T cells using a 2nd-generation packaging system. Briefly, cells were co-transfected with the shRNA lentiviral vector together with psPAX2 and pMD2.G at a plasmid ratio of 4:3:1, with a total of 8 µg DNA per 10 cm dish, using Lipofectamine 2000. After incubation at 37˚C for 6-8 h, the transfection medium was replaced with fresh complete medium. Lentiviral supernatants were collected 72 h post-transfection and filtered through a 0.45 µm membrane. ATDC5 cells were infected with the lentiviral supernatants at a multiplicity of infection of 10 in the presence of 10 µg/ml polybrene (Shanghai GenePharma Co., Ltd.). At 48 h, the medium was replaced, and downstream experiments were conducted 72 h after infection. Stably transduced cells were selected and maintained in the presence of 2 µg/ml puromycin.

Reverse transcription-quantitative PCR (RT-qPCR) analysis

Total RNA from transfected ATDC5 cells was isolated using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.) following the manufacturer's protocol. The mRNA was reverse transcribed into cDNA using the RevertAid RT Kit (Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. The target genes were amplified with an initial denaturation step at 95˚C for 10 min, followed by 40 amplification cycles consisting of 95˚C for 15 sec and 60˚C for 1 min. β-actin served as the internal reference for normalization. Relative mRNA expression levels were determined using the 2-ΔΔCq method as described previously (56,57). The primer sequences used were: LDLR forward, 5'-TCCAATCAATTCAGCTGTGG-3'; LDLR reverse, 5'-GAGCCATCTAGGCAATCTCG-3'; GAPDH forward, 5'-AACGACCCCTTCATTGAC-3'; and GAPDH reverse 5'-TCCACGACATACTCAGCAC-3'.

Western blotting analysis

ATDC5 chondrocytes were lysed in RIPA buffer (Beyotime Biotechnology) and the lysates were quantified using the BCA method. Equal amounts of proteins (30 µg) were subjected to 10% of sodium dodecyl sulfate-polyacrylamide gel electrophoresis and PVDF transfer. The membranes were blocked at room temperature in 5% milk for 1 h and incubated with specific primary antibodies at 4˚C overnight, including anti-LDLR (1:1,000; cat. no. ab30532; Abcam), anti-MYC (1:1,500; cat. no. ab32072; Abcam), anti-ADM (1:1,500; cat. no. ab190819; Abcam), anti-NFKBIA (1:2,000; cat. no. ab7217; Abcam), anti-MMP13 (1:1,000; cat. no. ab39012; Abcam) and anti-GAPDH (1:5,000; cat. no. ab8245; Abcam) antibodies. The membranes were subsequently washed with TBST (TBS containing 0.1% Tween-20) and incubated with HRP-conjugated anti-mouse (1:5,000; cat. no. #7076; CST Biological Reagents Co., Ltd.) or anti-rabbit (1;5,000; cat. no #7074; CST Biological Reagents Co., Ltd.) secondary antibodies at room temperature for 1 h. The protein bands were detected using enhanced chemiluminescence substrate (Shanghai Yeasen Biotechnology Co., Ltd.) (58).

Cell Counting Kit-8 (CCK-8) assays

Proliferation of transfected ATDC5 chondrocytes was assessed using a CCK-8 assay. Briefly, 100 µl cell suspension (4x103 cells/well) was seeded into 96-well plates and incubated at 37˚C for 24, 48 and 72 h. Subsequently, 10 µl CCK-8 (Shanghai Yeasen Biotechnology Co., Ltd.) reagent was added to each well and incubated for 2 h at 37˚C, after which the absorbance at 450 nm was measured (59).

Apoptosis analysis

After centrifugation at 150 x g for 5 min at 4˚C, ATDC5 cells were resuspended in 500 µl of binding buffer [Multi Sciences (Lianke) Biotech, Co., Ltd.] at 1x106 cells/ml. The cells were then stained with FITC-annexin V and propidium iodide (PI) by apoptosis detection kit [Multi Sciences (Lianke) Biotech, Co., Ltd.] according to the manufacturer's protocol and analyzed using FACSCalibur flow cytometry (BD Biosciences). Data were analyzed using FlowJo software (version 10.8.1; FlowJo LLC; BD Biosciences). Annexin V/PI quadrants were used to classify and quantify the viable, early apoptotic and late apoptotic/necrotic populations.

Statistical analysis

A total of six OA datasets were integrated and batch effects were removed to generate the final dataset that was used for analysis. IRDGs were extracted by comparing OA and control samples. Consensus clustering analysis was carried out on the basis of IRDG, followed by comparisons of gene expression, metabolic pathways, immune characteristics and drug sensitivity across different clusters. Combined LASSO regression analysis (60) and four ML algorithms were used to screen disease-associated biomarker genes, which were validated in three transcriptomic datasets. Target gene annotation was carried out using sc-transcriptomic data. All statistical analyses were carried out using the R programming language (version 4.1.3; Posit Software, PBC) (61) and Perl programming language (Strawberry Perl; version 5.30.0.1; Microsoft Corporation) (62). Perl scripts were primarily used for data integration and extraction, whereas R scripts were used for systematic data analysis and visualization. Experiments were repeated three times. Our data were presented as the mean ± SD. Differences between two groups were assessed using paired Student's t-tests, while comparisons among multiple groups were evaluated using one-way ANOVA with Tukey post-test. P<0.05 was considered to indicate a statistically significant difference.

Results

DEGs are identified across OA dataset from gene expression omnibus

Integration of the six datasets resulted in the identification of 8,238 co-expressed genes. Among these genes, 537 DEGs were identified between the OA and control samples, consisting of 204 genes with upregulated expression levels and 333 genes with downregulated expression levels. The six most significantly expressed DEGs were KLF transcription factor 9, MAF BZIP transcription factor F, JUND, GABA type A receptor associated protein like 1, heterogeneous nuclear ribonucleoprotein A1 and proline rich nuclear receptor coactivator 1 (Fig. 1A). By intersecting the 537 DEGs with 200 inflammation-associated genes, 11 IRDGs were identified (Fig. 1B), namely ADM, BTG anti-proliferation factor 2 (BTG2), cyclin dependent kinase inhibitor 1A (CDKN1A), IL-1 receptor type 1 (IL1R1), KLF transcription factor 6 (KLF6), LDLR, MYC, nicotinamide phosphoribosyltransferase (NAMPT), NFKBIA, presenilin 1 (PSEN1) and selectin E (SELE) (Fig. 1C). From the perspective of human chromosomes, BTG2 and SELE are located on chromosome 1. IL1R1, CDKN1A, NAMPT, MYC, KLF6, ADM and LDLR correspond to chromosomes 2, 6, 7, 8, 10, 11 and 19, respectively. PSEN1 and NFKBIA are located on chromosome 14 (Fig. 1D). Using the 11 IRDGs to construct a PPI network, five hub genes, NAMPT, NFKBIA, MYC, KLF6 and CDKN1A, were identified as central nodes in the network (Fig. 1E).

Figure 1 Differential gene expression analysis across Gene Expression Omnibus OA datasets. (A) Volcano plot of DEGs between patients with OA and controls. (B) Intersection of DEGs and inflammation-related genes. (C) IRDGs in patients with OA. (D) Chromosomal localization of IRDGs on 23 human chromosomes. (E) Protein-protein interaction networks constructed for the 11 IRDGs. DEGs, differentially expressed genes; IRDG, inflammation-related differentially expressed genes; OA, osteoarthritis. ADM, adrenomedullin; BTG2, B-cell translocation gene 2; CDKN1A, cyclin dependent kinase inhibitor 1A; IL1R1, interleukin 1 receptor type 1; KLF6, kruppel-like factor 6; LDLR, low density lipoprotein receptor; MYC, MYC proto-oncogene, NAMPT, nicotinamide phosphoribosyltransferase; NFKBIA, NFKB inhibitor alpha; SELE, selectin E.

OA subtypes are identified based on differentially expressed inflammatory genes

Consistent clustering consensus matrix based on the expression of 11 IRDGs classified the OA samples into two distinct clusters (clusters C1 and C2; Fig. 2A). The clustering result using k=2 was determined to be optimal on the basis of the cumulative distribution function curve (Fig. 2B) and the δ area curve (Fig. 2C). Cluster C1 contained 25 samples and the remaining 35 were assigned to cluster C2. The expression levels of IRDGs in different clusters are shown in Fig. 2D. PSEN1 was highly expressed in cluster C2 but at low levels in cluster C1, whereas the other IRDGs showed upregulated expression in cluster C1 compared with cluster C2. PCA demonstrated the distinct classification of clusters C1 and C2 (Fig. 2E). GSVA revealed the top 10 upregulated (‘prostate cancer’, ‘NOD like receptor signaling pathway’, ‘small cell lung cancer’, ‘pathways in cancer’, ‘adipocytokine signaling pathway’, ‘cytosolic DNA sensing pathway’, ‘nitrogen metabolism’, ‘apoptosis’, ‘JAK STAT signaling pathway’, ‘insulin signaling pathway’) and downregulated pathways (‘vasopressin regulated water reabsorption’, ‘long term potentiation’, ‘asthma’, ‘leishmania infection’, ‘biosynthesis of unsaturated fatty acids’, ‘gap junction’, ‘alpha linolenic acid metabolism’, ‘intestinal immune network for IgA production’, ‘ubiquitin mediated proteolysis’, ‘type 1 diabetes mellitus’) in cluster C2 compared with those in cluster C1 (Fig. 2F).

Figure 2 OA subtype classification based on differentially expressed inflammatory genes. (A) Consensus clustering matrix heatmap defining two clusters (k=2) and their correlation areas. (B) Consensus CDF curve demonstrating optimal clustering. (C) δ-area curve of consensus clustering in OA, indicating the relative change in the area under the CDF curve for each category number k compared with k-1. The horizontal axis represents the category number k, whereas the vertical axis represents the relative change in the area under the CDF curve. (D) Heatmap of the expression profiles of 11 inflammation-related differentially expressed genes in the two clusters. (E) PC analysis illustrating distinct immune subtypes in clusters C1 and C2. (F) GSVA pathway analysis comparing the top upregulated (red) and downregulated (blue) pathways between two distinct clusters. OA, osteoarthritis; CDF, cumulative distribution function; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; PC, principal component. ADM, adrenomedullin; BTG2, B-cell translocation gene 2; CDKN1A, cyclin dependent kinase inhibitor 1A; IL1R1, interleukin 1 receptor type 1; JAK, Janus kinase; KLF6, kruppel-like factor 6; LDLR, low density lipoprotein receptor; MYC, MYC proto-oncogene, NAMPT, nicotinamide phosphoribosyltransferase; NFKBIA, NFKB inhibitor alpha; STAT, signal transducer and activator of transcription; SELE, selectin E.

Different OA subtypes exhibit distinct immune characteristics

Only the stromal score differed among the immune, stromal and estimated scores between the two clusters (P=0.0011; Fig. 3A). Correlation analysis of the 22 immune cell types revealed primarily negative correlations, with neutrophils, activated mast cells and M2 macrophages exhibiting the strongest negative correlations with the other immune cell types (Fig. 3B). By contrast, memory B cells displayed the strongest positive correlation with the other immune cell types. Differential analysis of immune cell composition revealed that 10 immune cell types were differentially expressed (P<0.05) between clusters C1 and C2, with 8 exhibiting highly significant differences (P<0.001; Fig. 3C). In terms of immune functions, 7 immune functions were differentially enriched between the two clusters (P<0.05), with the most significantly altered immune functions being ‘Check-point’, ‘T cell co-inhibition’, chemokine receptor (‘CCR’) and antigen presenting cell ‘(APC) co-stimulation’ (Fig. 3D). The expression levels of immune checkpoints CD200 and CD160 significantly differed between the two clusters. Compared with cluster C1, cluster C2 exhibited increased CD160 expression but decreased CD200 expression (Fig. 3E).

Figure 3 Immune characteristics of different OA subtypes. (A) Defining two clusters (k=2) matrix score boxplots. (B) Correlation matrix of 22 immune cell types, stromal scores and immune scores in OA. Red indicates a positive correlation and blue indicates a negative correlation. (C) Differentially expressed immune cells between clusters C1 and C2 in OA. (D) Differentially expressed immune functions in clusters C1 and C2. (E) Expression levels of immune checkpoints in clusters C1 and C2. *P<0.05, **P<0.01 and ***P<0.001. OA, osteoarthritis; ns, not significant; NK, natural killer; HLA, Human leukocyte antigen; CCR, CC chemokine receptor; MHC, major histocompatibility complex; APC, antigen-presenting cell.

A ML model is constructed to identify disease signature genes

Before the ML model was constructed, the 11 IRDGs were screened using LASSO regression analysis (Fig. 4A), which narrowed the selection to eight disease signature genes, namely ADM, KLF6, LDLR, MYC, NAMPT, NFKBIA, PSEN1 and SELE (Fig. 4B). The residual box-line plots indicated that the SVM model had the smallest root-mean-square residuals (Fig. 4C). The reverse cumulative residual distribution likewise favored the SVM model (Fig. 4D). ROC curve analysis revealed that the SVM model achieved 98.5% accuracy, which demonstrated that it was the best-performing model for disease signature gene prediction (Fig. 4E). Calibration curves revealed that the predicted values were highly consistent with the actual results, although a high fluctuation at a 90% prediction probability was noted, requiring further investigation (Fig. 4F). A nomogram was constructed based on four key disease signature genes (MYC, ADM, LDLR and NFKBIA), which were selected by the SVM model, with each gene assigned a corresponding score to assess the risk of OA in individuals (Fig. 4G).

Figure 4 Machine-learning model construction. (A) LASSO regression analysis for preliminary screening of OA signature genes. (B) Selection of disease-related signature genes through LASSO regression. (C) Residual box line diagram for the GLM, RF, SVM and XGB models. The red dot indicates the root mean square of residuals. (D) Reverse cumulative distribution of residuals for different machine learning models. (E) Receiver operating curve analysis verifying model accuracy. (F) Calibration curve for assessing the accuracy of the nomogram. (G) Nomogram predicting OA incidence based on four signature genes. LASSO, least absolute shrinkage and selection operator; OA, osteoarthritis; GLM, generalized linear models; RF, random forests; SVM, support vector machines; XGB, extreme gradient boosting; ADM, adrenomedullin; LDLR, low-density lipoprotein receptor; NFKBIA, NF-κB inhibitor-α; KLF6, Krüppel-like factor 6; NAMPT, nicotinamide phosphoribosyltransferase; PSEN1, presenilin-1; SELE, selectin E.

Model accuracy is validated

The classification performance of the ML model was evaluated in three validation datasets based on the four identified genes (MYC, ADM, LDLR, and NFKBIA). ADM had a 100% classification accuracy in the GSE169077 dataset (Fig. 5A-D). Compared with the validation datasets, LDLR had the highest classification accuracy (72.7%) in the training set (Fig. 5E-H). MYC had an 83.3% classification accuracy in the GSE169077 dataset (Fig. 5I-L). NFKBIA had an 83.8% accuracy in the training set and an 83.3% accuracy in the GSE169077 validation dataset (Fig. 5M-P). Collectively, the classification accuracy for all the genes in the training set exceeded 70%. Among the test datasets, GSE169077 demonstrated the best classification performance, followed by GSE117999, whereas GSE82107 yielded the lowest accuracy. The differences in performance across datasets indicate that the classification ability of these genes is not completely uniform across cohorts. However, because the present study relied on retrospective public datasets with unavoidable differences in sample composition, sequencing platform, preprocessing workflow, and underlying biological heterogeneity, we could not determine the exact reason for the observed variation. Thus, our findings suggest potential predictive value of these genes, but further validation in larger, prospectively collected cohorts is required.

Figure 5 ROC curves of the four osteoarthritis disease signature genes in the training and test sets. ROC curves of ADM in the (A) training, (B) GSE82107, (C) GSE117999 and (D) GSE169077 datasets. ROC curves of LDLR in the (E) training, (F) GSE82107, (G) GSE117999 and (H) GSE169077 datasets. ROC curves of MYC in the (I) training, (J) GSE82107, (K) GSE117999 and (L) GSE169077 datasets. ROC curves of NFKBIA in the (M) training, (N) GSE82107, (O) GSE117999 and (P) GSE169077 datasets. ROC, receiver operating curve; ADM, adrenomedullin; LDLR, low-density lipoprotein receptor; NFKBIA, NF-κB inhibitor-α; AUC, area under the curve.

Four osteoarthritis signature genes are identified

OA sc sequencing datasets contained data on 1,464 chondrocytes derived from 10 patients with OA who underwent knee replacement surgery, with a total of 34 samples. Based on the intercellular coefficient of variation, the top 3,000 genes with the greatest variation among the 15,880 genes were extracted and included in the subsequent analysis (Fig. S1A-F). Sc data were analyzed using PCA for dimensionality reduction and t-SNE analysis for clustering (Fig. S2A-L). The expression patterns of four OA disease signature genes (ADM, LDLR, MYC and NFKBIA) were identified by the SVM model and analyzed using a GSE104782 sc dataset. The expression levels in the nine chondrocyte clusters revealed that NFKBIA had the highest expression levels, followed by MYC, ADM and LDLR (Fig. 6A). This finding was supported by t-SNE scatter plots (Fig. 6B-E). Cartilage synovial cell expression analysis revealed that when comparing the expression levels of ADM, LDLR, MYC and NFKBIA, NFKBIA was the most highly expressed and widely distributed, whereas LDLR was expressed at the lowest level (Fig. 6F-J).

Figure 6 Expression profile analysis of four osteoarthritis signature genes in the single cell RNA-sequencing data. (A) Bubble map displaying the expression levels of ADM, LDLR, MYC and NFKBIA genes in the nine clusters. A darker blue color indicates a higher expression level. t-SNE scatter plots illustrating the expression of (B) ADM, (C) LDLR, (D) MYC and (E) NFKBIA in chondrocytes. Red represents upregulated genes and green indicates downregulated genes. (F) Bubble map visualizing the expression levels of four genes in chondrocytes. Violin diagram showing the differential expression of (G) ADM, (H) LDLR, (I) MYC and (J) NFKBIA in chondrocytes. ADM, adrenomedullin; LDLR, low-density lipoprotein receptor; NFKBIA, NF-κB inhibitor-α; t-SNE, t-distributed stochastic neighbor embedding.

LDLR expression suppression inhibits cell viability and induces apoptosis

To investigate the function of LDLR in ATDC5 chondrocytes, the expression of LDLR was silenced using shRNA and overexpressed using LDLR cDNA. Knockdown using shRNA efficiently reduced LDLR expression levels at both the mRNA and protein levels compared with the control. Additionally, compared with the control, cDNA lentiviral delivery increased the expression levels of LDLR (Fig. 7A and B). Functionally, CCK-8 assays revealed that LDLR depletion decreased cell viability compared with the control. However, LDLR overexpression increased cell viability compared with the control (Fig. 7C). Furthermore, annexin V/PI staining revealed increased apoptosis after LDLR knockdown and reduced apoptosis after LDLR overexpression compared with the relevant controls (Fig. 7D and E). To investigate whether ADM, MYC and NFKBIA expression levels affected cell viability and apoptosis, their expression in ATDC5 cells was knocked down (Fig. S3A). Silencing ADM, MYC or NFKBIA reduced cell viability compared with the controls (Fig. S3B) and induced apoptosis in ATDC5 cells (Fig. S4). Collectively, these findings indicated that LDLR, ADM, MYC and NFKBIA may contribute to the regulation of ATDC5 cell viability and survival.

Figure 7 LDLR increases cell viability and suppresses apoptosis in ATDC5 chondrocytes. (A) Reverse transcription-quantitative PCR was carried out to measure the mRNA expression of LDLR in ATDC5 cells after LDLR modulation. (B) Western blotting was conducted to determine protein levels of LDLR in ATDC5 cells after LDLR modulation. (C) A Cell-Counting Kit 8 assay was carried out to evaluate ATDC5 cell viability after LDLR modulation. (D) Annexin V/PI flow cytometry was used to measure apoptosis in ATDC5 cells after LDLR modulation. (E) Quantitative results for apoptosis. (F) Western blotting was conducted to measure protein expression of MMP13 in ATDC5 cells after LDLR modulation. **P<0.01 and ***P<0.001. LDLR, low-density lipoprotein receptor; sh, short hairpin; EV, empty vector; OE, overexpression; OD, optical density.

LDLR suppression increases the expression of MMP13 in ATDC5 cells

Loss of LDLR function may cause a decreased LDL uptake by chondrocytes and an increased cholesterol deposition, which may activate inflammatory pathways (63). Excess cholesterol can activate the NF-κB pathway, induce the expression of inflammatory cytokines (such as IL-1β and TNF-α) and upregulate MMP13 expression levels, promoting the degradation of the cartilage matrix (64). To determine whether LDLR regulates the expression of MMP13 in ATDC5 cells, MMP13 expression levels were quantified following the manipulation of LDLR expression levels. Knockdown of LDLR expression levels through shLDLR transfection increased MMP13 expression levels compared with the control. However, overexpression of LDLR reduced MMP13 expression levels compared with the control (Fig. 7F). These data suggested that LDLR may limit cartilage matrix degradation by suppressing the expression of MMP13.

Discussion

In recent years, numerous studies have investigated OA pathogenesis using integrative bulk and scRNA-seq techniques. For example, a study by Peng et al (65) identifies 58 oxidative stress-related DEGs and seven diagnostic genes [stanniocalcin 2, lymphocyte specific protein 1, collagen type VI α-1 chain, Fos proto-oncogene, selenoprotein N, tumor protein p53 and heat shock protein family A (Hsp70) member 8] in OA through integration of GEO bulk and scRNA-seq. Furthermore, by integrating bulk and sc-transcriptomics with ML, a study by Sun et al (66) identifies the dysregulated expression of the hub genes S100 calcium binding protein (S100)-A9, phorbol-12-myristate-13-acetate-induced protein 1, ectodysplasin A2 receptor and fatty acid synthase, revealing that programmed cell death is associated with OA. Additionally, another study identifies four m7G regulators, eukaryotic translation initiation factor 1 (EIF1), JUND, nudix hydrolase 16-like 1 and nuclear cap binding protein subunit 1, as diagnostic biomarkers for OA (67).

EIF1/JUND suppression induces MMP13 upregulation, collagen II downregulation and OA progression (67). An additional study mapped endoplasmic-reticulum-stress heterogeneity in osteoarthritic chondrocytes and identified insulin like growth factor binding protein 3 and S100A4 as key diagnostic regulators (68). Tryptophan metabolism, through the activation of the indoleamine 2,3-dioxygenase 1/tryptophan 2,3-dioxygenase-induced kynurenine pathway and the accumulation of 3-hydroxykynurenine/quinolinic acid, promotes cartilage degradation in OA chondrocytes (69). In the present study, 11 IRDGs were identified in the OA and normal samples. The regulatory network, composed of five hub genes (NAMPT, NFKBIA, MYC, KLF6 and CDKN1A), was highly coincident with the NF-κB/STAT3 inflammatory signaling axis and the specific expression of PSEN1 in cluster C2 suggested that it may exacerbate cartilage degradation through mitochondrial dysfunction. Notably, PSEN1 has been studied in association with Alzheimer's disease (AD) (70-72). Mutations in PSEN1 reduce barrier function, impair glucose metabolism and decrease drug efflux pump activity (73). Given that both OA and AD are degenerative diseases with a high prevalence in elderly individuals, future studies may benefit from investigating whether PSEN1 dysregulation contributes to OA pathogenesis.

In the inflammatory milieu, immune cell types were predominantly negatively correlated with one another. Significant differences were observed in numerous immune cell populations (including memory B cells, T regulatory cells, activated natural killer cells, M0 macrophages, resting dendritic cells, resting mast cells, activated mast cells and neutrophils), four immune functions (checkpoint, T cell co-inhibition, CCR and APC co-stimulation) and one immune checkpoint (CD200). The present findings regarding immune cell infiltration and immune function align with those of a previous OA immune infiltration study (74). In patients with posttraumatic OA, expression of the cell surface glycoprotein CD200 receptor 1 is markedly increased and contributes to cartilage degeneration (75), suggesting that CD200 may be a key inflammation-related biomarker in patients with OA. Mechanistically, the Akt inhibitor MK-2206 has anti-autophagic, proapoptotic, prosenescent and procatabolic effects in human intervertebral discs. Conversely, Akt activation in intervertebral disc cells protects against inflammation through temsirolimus treatment (76). However, further research is needed to investigate the potential role of these drugs in OA treatment.

Using LASSO regression analysis in combination with four ML algorithms (GLM, RF, SVM and XGB), MYC, ADM, LDLR and NFKBIA were identified as inflammation-related biomarkers of OA. One study suggests that MYC mediates aberrant gene expression associated with OA and may serve as a potential therapeutic target (35). Bioinformatics of OA synovium shows elevated expression of MYC, JUN and VEGFA, which is associated with severity. Furthermore, MYC/VEGFA may induce p38-MAPK/JUN-mediated inflammation, proliferation and angiogenesis (77). Borojoa iridoid glycoside enhanced chondrocyte proliferation and reduced apoptosis through MYC upregulation (78). In ATDC5 cells, the miR-195-5p inhibitor restores proliferation and reduces apoptosis and inflammation through the activation of Wnt/β-catenin and NF-κB pathways, including MYC, cyclin D1 and p-p65(79). ATDC5 cells lacking nuclear factor of activated T (NFATc)-1 and NFATc2 have an increased proliferation and an upregulated expression of MYC (80). In IL-1β-induced ATDC5 cells, overexpression of long non-coding RNA colorectal neoplasia differentially expressed enhances proliferation and reduces apoptosis, inflammation and extracellular matrix (ECM) degradation, while IL-1β activates Wnt/β-catenin signaling with increases in the expression of β-catenin, MYC and cyclin D1(81). MYC overexpression enhances chondrocyte viability and suppresses ferroptosis. In vivo, MYC overexpression mitigates cartilage damage in OA mouse models (82). These findings suggested that MYC may promote chondrocyte viability and inhibit apoptosis. Consistently, the present data showed that MYC silencing reduced viability and increased apoptosis in ATDC5 cells.

A number of studies suggest that ADM may serve as a useful diagnostic biomarker for OA (83-86). For example, microarray profiling of OA bone identifies 150 DEGs, including ADM, implicating Wnt/TGF-β/bone morphogenetic protein and osteoclast pathways in altered bone remodeling and sex differences (87). Shear stress changed the transcriptional profile of osteoarthritic chondrocytes, with microarray analysis identifying 14 upregulated and 6 downregulated genes, including ADM, that may increase cell survival and matrix homeostasis (88). A transcriptome-wide association study, integrating genome-wide association studies and expression profiles, identifies shared OA genes and pathways, including ADM, FSTL1 genes and the PI3K/AKT pathway. ADM is a notable candidate, highlighting PI3K/AKT and immune-associated mechanisms, such as inflammatory signaling, antigen presentation, chemokine activity, and innate immune responses (89). Through integrating DEGs and weighted gene co-expression network analysis, a previous study identifies nuclear factor IL-3 regulated, ADM and osteoglycin as early OA biomarkers, associating MYC targets and NF-κB signaling with inflammatory pathways, including ROS and TNF-a (90). Using GSVA/single sample GSEA, hypoxia is implicated in OA and immune infiltration. A seven-gene index including ADM achieved high diagnostic accuracy (91). Consistently, the present study identified that ADM may be an inflammation-related OA biomarker. ADM prevents fibroblast-like synovial apoptosis in patients with RA by upregulating calcitonin receptor like receptor/receptor activity modifying protein 2 expression (92). In line with this finding, the present results showed that silencing ADM reduced cell viability and induced apoptosis in ATDC5 cells.

NFKBIA, an immune-related gene, has high predictive performance in OA diagnostic models (93-95). In the Han Chinese, SNPs in TLR10 and NFKBIA are associated with hip OA risk (96). Microarray analysis of joint fibroblast-like synoviocytes identifies 276 OA DEGs. The OA-specific hubs genes SELE, SERPINE1 and NFKBIA implicated that TNF signaling may contribute to OA inflammation by mediating the pro-inflammatory activity of fibroblast-like synoviocytes in the synovium (97).

Safflower yellow, a natural pigment extracted from safflower, protects chondrocytes and suppresses OA inflammation by inhibiting TNF-α-induced NF-κB activation, reducing MMP-13 expression levels and improving cartilage degeneration (98). Integrating OA synovium DEGs with aging genes, ML identifies activating transcription factor 3, KLF transcription factor 4, NFKBIA and superoxide dismutase 2 as downregulated diagnostic biomarkers associated with immune infiltration (99). In addition, tRNA-derived fragment 16 is upregulated in OA and worsens disease by targeting AlkB homolog 5 RNA demethylase, reducing NFKBIA mRNA stability, activating NF-κB and promoting inflammation and ECM degradation (100). Microarray and CIBERSORT analyses identifies macrophages as dominant infiltrates in OA, defines two immune-related OA subgroups and highlights hub genes including NFKBIA and MYC (101). The present study identified NFKBIA as a potential OA biomarker.

LDLR serves a pivotal role in OA pathogenesis. In a murine OA model, LDL accumulation through a cholesterol-rich diet or LDLR deficiency promotes synovial macrophage oxLDL uptake, increases S100A8 and TGF-β signaling and markedly increases osteophyte and enthesophyte formation (102). Through the use of drug-target Mendelian randomization, genetically proxied LDL reduction through 3-hydroxy-3-methylglutaryl-CoA reductase, proprotein convertase subtilisin/kexin type 9 and LDLR is associated with a reduced risk of OA (103). By integrating multi-tissue transcriptomics with ML and scRNA-seq, a previous study identifies five peripheral blood biomarkers, including LDLR, that may be used to predict OA, reveal age-specific expression profiles and B cell remodeling, as well as associations with immune infiltration and inflammatory cytokines (104). In the present study, it was revealed that LDLR increased cell viability, inhibited apoptosis and repressed the expression of MMP13 in ATDC5 cells, indicating a potential role for LDLR in OA pathogenesis. Intracellular cholesterol accumulation may be one of the upstream events involved in activating inflammatory pathways such as NF-κB, which is a key transcription factor driving the expression of matrix degrading enzymes such as MMP13. The results of the present study supported the concept that there may be a functional association between lipid metabolism disorders and the joint inflammation catabolic axis.

In conclusion, the present study integrated six GEO bulk transcriptomic cohorts with single-cell chondrocyte profiling and multi-algorithm ML to define inflammation-driven OA subtypes. A total of four key biomarkers were identified (ADM, LDLR, MYC and NFKBIA) and their expression was validated using scRNA-seq data. In addition, the present study added functional perturbation of ADM, LDLR, MYC and NFKBIA in ATDC5 cells. However, a number of limitations should be acknowledged. Firstly, the proposed mechanistic axis, such as the cholesterol/NF-κB/MMP13 axis, was primarily hypothesis-driven and was not directly investigated with causal experiments (such as pathway perturbation or rescue assays). Secondly, the biomarkers investigated in the present study relied on computational integration of public transcriptomic datasets, which may have introduced platform- and batch-related biases, residual confounding and model overfitting. Thirdly, the sample size was modest, which may limit the generalizability across diverse OA phenotypes and clinical stages. Furthermore, functional validation was carried out in ATDC5 cells, a mouse chondrogenic cell line, and potential differences compared with human OA chondrocytes and joint tissues were not evaluated. Subsequently, future work should prioritize validation of the present findings using multicenter human cohorts, including synovium and cartilage samples, longitudinal follow-up and treatment-response data. Mechanistic experiments using human primary chondrocytes, organoid models or mouse models with targeted perturbation of cholesterol handling and NF-κB activity should be carried out to determine causality. Lastly, the SVM-based model and nomogram may support earlier OA detection and immune-inflammation stratification from synovial or cartilage expression data, enabling risk scoring and cohort enrichment for trials. The LDLR and MMP13 could be targets and potential biomarkers for OA treatment and monitor response.

Supplementary Material

Quality control and filtering of single cell RNA-sequencing data. (A) Number of genes detected in 34 osteoarthritis samples. (B) Sequencing depth across all samples. (C) Mitochondria gene content in each sample. High levels meant low cell activity. (D) Association between sequencing depth and mitochondrial gene content. (E) A positive correlation between sequencing depth and total intracellular sequences (R=0.25). (F) Volcano plot illustrating top 3,000 highly variable genes. Red indicates high variation. The names of the top 10 genes are presented.

PCA downscaling and t-SNE clustering of scRNA-seq data. Bubble plot showing the top 20 marker genes for each of the PCA clusters (A) PC 1, (B) 2, (C) 3 and (D) 4. Heatmaps of marker genes in each of the PCA clusters (E) PC 1, (F) 2, (G) 3 and (H) 4. Yellow represents upregulated genes, while purple indicates downregulated genes. The top 30 genes are labeled. (I) PCA projection of scRNA-seq data for preliminary dimensionality reduction. (J) P-values for PCA components with the top 15 PCs selected for subsequent analysis. (K) t-SNE was used for clustering of osteoarthritis cells. (L) Cell type annotation of clusters. Nine chondrocyte clusters were determined. PC, principal component; PCA, PC analysis; t-SNE, t-distributed stochastic neighbor embedding; scRNA-seq, single cell RNA-sequencing.

Silencing MYC, ADM or NFKBIA reduces cell viability in ATDC5 chondrocytes. (A) Western blotting was carried out to detect the protein levels of MYC, ADM and NFKBIA in ATDC5 cells after shRNA transfection. (B) A Cell Counting Kit-8 assay was carried out to evaluate cell viability in ATDC5 cells after shRNA transfection. ***P<0.001. ADM, adrenomedullin; NFKBIA, NF-κB inhibitor-α; OD, optical density; sh, short hairpin.

Silencing MYC, ADM and NFKBIA induces apoptosis in ATDC5 chondrocytes. (A) Annexin V/PI flow cytometry was used to measure the apoptosis of ATDC5 cells after shRNA transfection. (B) Quantitative results for apoptosis. ***P<0.001. ADM, adrenomedullin; NFKBIA, NF-κB inhibitor-α; sh, short hairpin.

Patient clinical information.

Acknowledgements

Not applicable.

Funding

Funding: The present study was supported by The Jiande City Science and Technology Development Project (grant nos. 2023SJZX13 and 2024SJZX08), The Zhejiang Provincial Medical and Health Science and Technology Project (grant no. 2022RC244) and The Hangzhou Agricultural and Social Development Project (grant no. 20241029Y206).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

JX performed experiments, analyzed the data and wrote the manuscript. WQW, XTL, SWX and BZ performed experiments and analyzed the data. XL conceived and designed the study, drafted the manuscript and supervised the study. JX and XL confirm the authenticity of all the raw data. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References