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.

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 log₂-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.

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.

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.

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).

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% CO₂ 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.

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'.

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).

Proliferation of transfected ATDC5 chondrocytes was assessed using a CCK-8 assay. Briefly, 100 µl cell suspension (4x10³ 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).

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 1x10⁶ 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.

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.

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).

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).

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).

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).

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.

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).

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.

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.