Gene expression profiles were obtained from the Gene Expression Omnibus (GEO) database. Five publicly available microarray datasets were included: GSE55235, GSE55457, and GSE82107 (derived from synovial tissue), and GSE51588 and GSE169077 (derived from articular cartilage). In total, these datasets comprised 76 osteoarthritis (OA) samples and 42 healthy control samples (As shown in Table 1). A list of 1,584 coagulation-related genes (CRGs) was retrieved from the GeneCards database (https://www.genecards.org/), including genes with a “Relevance score” > 1 that encode proteins involved in coagulation pathways. This threshold was applied to select genes with at least one or more lines of established evidence supporting their involvement in coagulation pathways, thereby constructing a biologically relevant and high-confidence candidate gene set for subsequent analysis (22, 23). To construct analysis cohorts that are representative of the joint microenvironment and to control for potential bias arising from tissue heterogeneity, we designed both the training set and the independent test set to contain a mixture of cartilage and synovial tissue samples. For cohort construction, expression matrices from GSE51588 (cartilage), GSE55457 (synovium), and GSE55235 (synovium) were merged and integrated using the removeBatchEffect function from the limma R package to form the training set. Similarly, GSE169077 (cartilage) and GSE82107 (synovium) were merged and batch-corrected using the same procedure to form the independent test set. The successful removal of batch effects was confirmed by comparative box plots of normalized expression values across the integrated datasets.

Using the R statistical tool “limma,” we carried out differential expression analysis comparing healthy and normal control samples with OA samples in the training set to identify differentially expressed genes (DEGs). Volcano plot of DEGs was drawn. The CRGs and DEGs were taken to intersect to obtain coagulation-related differentially expressed genes (DECRGs). Protein interaction analysis and visualization of DECRGs were performed using string database and cytoscape software.

We employed three distinct machine learning algorithms—Generalized Linear Model (GLM), Support Vector Machine (SVM), and Random Forest (RF)—to screen differentially expressed coagulation-related genes (DECRGs). To identify the most robust signals, the top 50 genes ranked by importance from each algorithm were selected for intersection. This cutoff was chosen because a clear inflection point was observed around the top 50 genes in the importance ranking curves across all methods, beyond which scores plateaued. This approach captures the most discriminatory features while focusing the consensus analysis on high-potential candidates (24). The protein-protein interaction (PPI) network centrality measures (Betweenness, MNC, Eccentricity) were calculated, and genes ranking in the top 50% for each measure were selected to capture topologically important nodes. The final set of core diagnostic genes was determined by taking the intersection of the gene sets derived from the machine learning consensus and the PPI network centrality analysis. This stringent, multi-perspective approach ensured the identification of robust biomarkers.

Violin plots for differential expression were generated using ggstatsplot. The area under the receiver operating characteristic curve (AUC) was calculated for each candidate gene using the ROCR package. Genes with an AUC greater than 0.7 were selected to build the diagnostic model. This widely adopted threshold ensures that only features with at least “acceptable” discriminatory power (AUC 0.7-0.8) are included, enhancing the model’s foundational validity (25–27). A multivariate logistic regression model was then built using the glmnet package with these selected genes. The model’s performance was evaluated by computing its AUC in both training and test sets. Finally, a nomogram was developed and calibrated using the rmda package, and decision curve analysis was performed to confirm its clinical utility.

To investigate the potential of coagulation-related genes in stratifying osteoarthritis (OA), we performed unsupervised consensus clustering on OA samples. In the training set, clustering was conducted using the expression profiles of the two core genes (COL3A1 and MMP1) with the ConsensusClusterPlus R package, evaluating cluster numbers (K) from 2 to 9 to define preliminary subtypes. To assess the robustness of this classification approach, the same consensus clustering procedure (setting K = 2-9) was independently applied to the OA samples in the test set using the identical gene features. The resultant subtype labels from both cohorts were subsequently used as grouping variables for downstream comparative analyses, including immune microenvironment characterization.

We identified differentially expressed genes (DEGs) for the new OA classification and performed enrichment analyses on these DEGs, including GSEA, KEGG, and GO analyses. A corrected p-value < 0.05 was used as the filtering criterion. This approach further elucidated the pathological mechanisms underlying stratified management of OA. We also performed subclustering based on the expression levels of the most critical core genes, employing the same methodology to explore the pathological mechanisms by which genes influence disease progression.

To comprehensively characterize the immune microenvironment across different OA subtypes, our analysis encompassed immune cell infiltration profiling and immune checkpoint gene expression assessment.

For primary cell isolation, articular cartilage from 5-day-old wild-type mice was used in this study. All procedures were performed under sterile conditions in a biosafety cabinet. Using ophthalmic scissors, the separated cartilage was carefully cut into tiny tissue pieces on a sterile Petri plate. The fragments were pretreated for tissue by digesting them with a 0.25% trypsin solution while being constantly stirred(30 minutes, 37°C). This was followed by further enzymatic digestion in DMEM over six hours. Following two rounds of centrifugation to purify the isolated cells, the pellet was resuspended in full DMEM/F12 media containing 1% P/S with 10% FBS. After that, cell suspension was seeded into standard T25 flasks and kept in an incubator with humidity for primary growth (5% CO2, 37°C), changing the medium every 48 hours. To ensure functional stability, only passage 1 to passage 2 (P1-P2) cells were used in subsequent experiments. During subculturing, cell adherence and morphological traits were regularly observed.

The mouse COL3A1 overexpression vector (pCMV-COL3A1-Neo, pCOL3A1) and its negative control (pNC) were constructed by Tsingke Biotechnology. Plasmids were transfected into chondrocytes using LipoLeagene2000 according to the manufacturer’s instructions, with untreated cells serving as a blank control. To confirm successful overexpression, qRT-PCR was performed on cells from each experimental group. Total RNA was isolated under equal input conditions using the High-Purity RNA Isolation Kit and reverse-transcribed into cDNA using the High-Capacity cDNA Reverse Transcription Kit. The following primer sets were used for qRT-PCR:

To mimic an osteoarthritic microenvironment, chondrocytes from both the empty vector-transfected group (pNC) and the COL3A1-overexpressing group (pCOL3A1) were given 10 ng/mL of IL-1β for 12 hours after transfection, designated as pNC+IL-1β and pCOL3A1+IL-1β, respectively. These groups, along with the untreated pNC group, were subsequently subjected to further cellular assays. After equally seeding 6-well plates with the chondrocyte suspension, the plates were incubated for a full day at 37°C with 5% CO2. Following the removal of growth medium, we fixed the cells for ten minutes at room temperature using 4% paraformaldehyde (PFA), and then we washed them three times with PBS to get rid of any remaining fixative. In accordance with the EdU assay kit’s manufacturer’s instructions, cells were stained with DAPI solution for nuclear visualization after being incubated with the EdU reaction mixture while shielded from light. Finally, after washing the samples with PBS, we observed them at the fluorescence microscope. The assay was performed in three independent biological replicates, with three technical replicates per experiment. EdU-positive cells were counted in three randomly selected fields per well using ImageJ software. Data are presented as mean ± SD.

Following their respective protocols, chondrocytes from each experimental group were seeded into 12-well plates and maintained in culture for 24 hours. To prepare for staining, the existing medium was removed, after which the cells in each well were exposed to pre-warmed JC-1 working solution (5 μM final concentration) and kept at 37°C in the absence of light for a 30-minute incubation period. Following the staining procedure, cells underwent three gentle washes with ice-cold PBS to eliminate any unbound dye. Visualization was carried out via an inverted fluorescence microscope configured for dual-channel acquisition: the red fluorescence from JC-1 aggregates was recorded with a 590 nm emission filter, and the green signal from monomers was captured with a 530 nm filter. The ratio of red-to-green fluorescence intensity was quantified using ImageJ software. Data were derived from three independent biological replicates.

To precisely identify cell surface chondroitin sulfate for Toluidine Blue staining, cells were treated with a 0.1% working solution over sixty minutes in a dark-controlled condition. Safranin O staining (0.02% working solution, 45 minutes) was employed to evaluate glycosaminoglycan distribution. Alcian Blue staining utilized a 1% solution (50-minute incubation) to visualize acidic mucopolysaccharide components. Following each staining step, we washed samples three times with PBS under gentle agitation (5 minutes per wash) to eliminate residual dye. Finally, we dehydrated specimens through a graded ethanol series, and we used inverted microscope to get broad-field pictures. Quantitative analysis of all images was performed under blinded conditions. The mean optical density (OD) was measured from three random fields per sample using ImageJ software. Data are expressed as the mean ± SD of at least three independent biological replicates (N = 3).

The DCFH-DA probe (Servicebio) was used to label cells in order to determine total intracellular ROS. The cell suspension was combined with the probe’s working solution (1:1000 dilution) and cultured over 30 minutes at 37°C in a light-protected environment. The unbound probe was then eliminated by rinsing it three times in HBSS buffer that had been heated beforehand. Fluorescence was visualized with the fluorescence microscope. Quantitative image analysis was performed by calculating the mean fluorescence intensity (MFI) with ImageJ software, with at least three biological replicates per treatment group included for statistical analysis. Quantitative image analysis was performed by calculating the mean fluorescence intensity using ImageJ software. Each treatment group included at least three biological replicates, with three technical replicate samples analyzed per experiment.

In this study, we evaluated cellular senescence using β-galactosidase staining. Following 3 full days of treatment with βOHB and TBHP, a solution containing 2% glutaraldehyde and 0.2% formaldehyde was used to fix chondrocytes from various experimental groups over fifteen minutes at room temperature. After completing the above steps, cells were incubated with working solution (Beyotime) at 37°C for 16 hours under light-protected conditions. After staining, samples underwent three graded PBS washes. Senescence-associated β-galactosidase activity was visualized by detecting green precipitates. The percentage of SA-β-Gal-positive cells was counted in three randomly selected fields per sample using ImageJ software. The experiment was independently repeated three times (N = 3).

R (with 4.4.2 version) was used to conduct statistical analysis. All cell-based experiments included at least three biological replicates. We sourced the transcriptomic data from the Gene Expression Omnibus repository (GEO; https://www.ncbi.nlm.nih.gov/geo/). For comparisons across groups, a two-tailed Student’s t-test was applied, setting the significance level at P < 0.05. *, **, and *** indicate p-value <0.05, p-value <0.01, and p-value <0.001, respectively, and NS indicates no statistically significant difference. The sample size for each group was three, and the experiments for each sample were performed in three biological replicates.

Figures 1A, B display the expression box plots for the training set before or after removing batch effect, while Figures 1C, D show corresponding plots for the test set. We observed that batch effects were successfully removed in both datasets. Figure 1E presents the DEGs across different samples, genes with elevated expression in OA tissues are shown in red, while those more abundant in normal samples are depicted in blue. We identified 73 differentially expressed coagulation-related genes (DECRGs) by taking the intersection CRGs and DEGs, as shown in Figure 1F. A protein-protein interaction (PPI) network for these 73 DECRGs was constructed and is displayed in Figure 1G. The positive or negative correlation relationships among these DECRGs are illustrated in Figure 1H. To evaluate the disease-predictive potential of the 73 DECRGs, we employed three machine learning methodologies: RF, SVM, and GLM. The results of the feature importance ranking from these methods are shown in Figures 2A–C, respectively. Additionally, we calculated and ranked the DECRGs based on EcCentricity, Betweenness, and MNC values. By taking the intersection of the top 50% of genes from all six ranking methods, we identified five key genes crucial for OA: COL3A1, CXCL8, IL1B, MMP1, and MMP3, as summarized in Figure 2D. The correlation relationships among these five key genes are visualized in a correlation chord diagram (Figure 2E) and a node-link diagram (Figure 2F). To further explore the disease predictive value of these key genes, we generated violin plots comparing their expression levels of OA samples with healthy control samples. As displayed in Figures 2G–K, OA samples exhibited markedly elevated expression levels of COL3A1, MMP1, and MMP3 relative to healthy controls. We plotted ROC for evaluating the independent diagnostic accuracy of each of the five key genes for OA. The results are presented in Figures 2L through 2P. Both COL3A1 and MMP1 demonstrated discriminative ability, with AUC values exceeding 0.7, suggesting that COL3A1 and MMP1 hold superior diagnostic potential for OA among these evaluated gene candidates.

We utilized COL3A1 and MMP1, which had the highest AUC values, to perform multivariate logistic regression for constructing an OA diagnostic model. The model is represented as: Riskscore = (MMP1 expression × 0.3893437) + (COL3A1 expression × 0.7646114). Each sample was assigned a risk score estimating its probability of having OA. For the training cohort, a risk score-based ROC curve was developed to assess OA predictive accuracy. This model demonstrated strong discriminative power with an AUC of 0.813, validating its diagnostic reliability. Unsupervised clustering was applied to training set OA samples, utilizing expression patterns of MMP1 and COL3A1 for classification. Our results indicated optimal clustering performance when the cluster number (K) was set to 2, as shown in Figures 3A–C. Additionally, we built a nomogram for OA prediction according to COL3A1 and MMP1 expression (Figure 3E). The calibration curve (Figure 3F) and DCA (Figure 3G) of the nomogram both demonstrated its high predictive potential. To further confirm the stability of our established OA subtyping and diagnostic model, we evaluated it in the test set. ROC curves were plotted for the individual risk genes (COL3A1 and MMP1) and the combined diagnostic model, yielding AUC of 0.766, 0.703, 0.766 (Figures 4A–C). Furthermore, clustering based on COL3A1 and MMP1 expression in the test set also showed excellent performance, as presented in Figures 4D–G.

Differential expression analysis was conducted between C1 and C2 OA subtypes to elucidate their distinct pathological mechanisms and evaluate variations in biological functions and molecular pathways, we also performed further enrichment analysis on these DEGs. The enrichment results of DEGs in MF, CC, and BP are displayed in Figure 5A. These genes showed significant enrichment in pathways associated with cartilage degeneration including collagen catabolic process, extracellular matrix containing collagen, collagen trimer, fibrillar collagen trimer, glycosaminoglycan binding, peptide binding, and serine hydrolase activity. Additionally, GO analysis revealed immune-related pathways such as neutrophil chemotaxis, leukocyte migration regulation, and neutrophil migration. In Figure 5B, the GO findings are displayed as circle diagrams. KEGG enrichment results are shown in Figures 5C, D, and the genes were significantly enriched in protein digestion and absorption, Rheumatoid arthritis, IL-17 signaling pathway, Lipid and atherosclerosis, TNF, AGE-RAGE, Relaxin, PPAR, Chemokine signaling pathway. With reference to the GSEA enrichment analysis, on the other hand, indicated that these genes appeared substantially abundant in pathways including GNF2 CD48, adaptive immune response, cartilage development, chondrocyte differentiation, leukocyte chemotaxis, glycosaminoglycan binding, coagulation, osteoarthritis, collagen degradation, TNF signaling via NFKB, as shown in Figure 5E. These pathways are not only associated with coagulation and osteoarthritis development, but also contain immune infiltration-related pathways. As a result, we postulated that there are variations in immune infiltration among various subtypes, which merits further investigation as a treatment approach for osteoarthritis.

Enrichment results indicated potential immune pathway variations across OA subtypes. Subsequent CIBERSORT and immune infiltration analyses were therefore conducted to characterize these immunological differences. Immune infiltration profiling across all samples (Figure 6A) revealed macrophages as the predominant immune cell population. Subtype comparisons in Figure 6B demonstrated markedly distinct immune activity profiles, including NK cell infiltration, antigen-presenting cell co-stimulation, type II interferon responses, and immune checkpoint expression. Figure 6C revealed substantial disparities in immune checkpoint expression (CD48, CD44, CD86, CD28, TNFSF4, TNFRSF25) across OA subtypes. Figure 6D demonstrated distinct immune infiltration patterns between C1 and C2 subtypes, with C1 exhibiting moderately elevated immune infiltration levels and significantly higher expression of both COL3A1 and MMP1 compared to C2.

We next examined how pivotal genes correlate with immune infiltration patterns in OA. We selected COL3A1, which had a higher risk coefficient (0.7646114) and a higher AUC value (0.766) in the test set, for subsequent analysis. After applying Benjamini-Hochberg false discovery rate (FDR) correction for multiple testing, COL3A1 expression levels showed significant positive correlations with CD44, CD86, and Type II IFN response, and a significant negative correlation with Inflammation-promoting (FDR-adjusted P < 0.05; Figures 7A–D). These results suggest that COL3A1 is associated with a distinct immune-inflamed phenotype. This finding is consistent with our observation that COL3A1 expression was substantially elevated in the C1 OA subtype, which also demonstrated markedly higher levels of immune infiltration compared to the C2 subtype. Figure 7E demonstrated the GSEA enrichment results of COL3A1, which showed that COL3A1 was associated with pathways including GNF2 TNFRSF1B, chondrocyte development, endochondral bone morphogenesis, response to interferon gamma, coagulation, interferon gamma response, arthritis, osteoarthritis, collagen formation, and inflammatory response. These pathways are collectively associated with osteoarthritis progression, coagulation, and immunity, suggesting that COL3A1 may affect OA progression by regulating biological processes such as chondrocyte development, subchondral bone formation, and immune-metabolic crosstalk.

To investigate the impact of COL3A1 on the chondrocyte phenotype, we established a COL3A1-overexpressing chondrocyte model (pCOL3A1) and a control (pNC) via plasmid transfection, with successful overexpression confirmed by PCR (Figure 8A). Subsequently, both pNC and pCOL3A1 cells were treated with IL-1β to simulate the osteoarthritic inflammatory environment. Cytochemical staining with Safranin O, Toluidine Blue, and Alcian Blue demonstrated that IL-1β stimulation markedly reduced the content of glycosaminoglycans and proteoglycans, indicating impaired extracellular matrix synthesis (Figure 8B). Quantitative analysis confirmed that this IL-1β-induced impairment was significantly exacerbated by COL3A1 overexpression (Figures 8C–E). Furthermore, the EdU assay revealed that IL-1β inhibited chondrocyte proliferation, an effect that was further enhanced in COL3A1-overexpressing cells, as shown by representative images (Figure 8F) and quantitative analysis (Figure 8G). Collectively, these results indicate that COL3A1 overexpression suppresses both the anabolic function and proliferative capacity of chondrocytes under inflammatory OA conditions.

To further elucidate the mechanisms by which COL3A1 affects chondrocyte function, we assessed mitochondrial function, oxidative stress, and senescence-related markers. Analysis of mitochondrial membrane potential via JC-1 staining revealed that IL-1β treatment led to a marked decrease in potential, an effect that was significantly exacerbated by COL3A1 overexpression, as shown in representative images (Figure 9A) and confirmed by quantitative analysis (Figure 9B). Similarly, measurement of intracellular reactive oxygen species (ROS) levels demonstrated that IL-1β stimulation significantly elevated oxidative stress, and COL3A1 overexpression further amplified this ROS production (Figures 9C, D). Finally, assessment of cellular senescence via SA-β-gal staining confirmed that IL-1β promoted chondrocyte senescence, a pro-senescent effect that was markedly enhanced by COL3A1 overexpression (Figures 9E, F). In summary, these results demonstrate that under osteoarthritic conditions, COL3A1 overexpression exacerbates chondrocyte degeneration by impairing mitochondrial function, intensifying oxidative stress, and promoting cellular senescence.