Comprehensive Analysis of Programmed Cell Death-Related Genes in Diagnosis and Synovitis During Osteoarthritis Development: Based on Bulk and Single-Cell RNA Sequencing Data

Introduction

Osteoarthritis (OA) is a common degenerative joint disease characterized by cartilage degradation accompanied by changes in joint structure such as synovitis, subchondral bone sclerosis, osteophyte formation, meniscus tears, and infrapatellar fat pad inflammation. Risk factors for OA include aging, female sex, obesity, genetic predisposition, and previous joint injuries, such as tibial plateau fractures, knee ligament injuries, and meniscal tears. According to studies, the prevalence of OA among adults in the worldwide reach 500 million people, which imposes a heavy burden on the healthcare system, society, and the economy.^1,2 Synovium serves as a bridge for cartilage nutrient supply and waste removal.³ The anatomical functional unit formed by synovium and IPFP plays an important role in the stability of knee joint structure, the regulation of inflammation and cartilage homeostasis.⁴ The chronic inflammatory pannus formed by synovial lesions in OA can lead to the degradation of cartilage matrix, and collagen fragments along with proteoglycans are released into the synovial fluid. This process triggers an inflammatory response in the synovial tissue, which further accelerates cartilage degeneration, narrowing of the joint space.^5,6 Therefore, clarifying the molecular mechanism of OA synovitis is crucial for OA early diagnosis, treatment, and prevention of further joint damage.

Programmed cell death (PCD) is a sophisticated, genetically encoded process that governs cellular self-termination.⁷ Prior research has demonstrated a strong connection between programmed cell death, such as apoptosis, autophagy, necroptosis, pyroptosis, ferroptosis, and cuproptosis, and the development of OA.⁸ Apoptosis, a highly regulated and active form of cell death, has been identified as a pivotal factor in the cartilage degeneration characteristic of OA.^9,10 Autophagy, a lysosomal degradation process, operates protectively within chondrocytes, preserving cellular integrity by clearing misfolded proteins and damaged organelles.^11,12 Necroptosis, regulated by receptor-interacting protein kinases RIPK1 and RIPK3, triggers the death of synovial fibroblasts, contributing to the development of synovitis in OA.^13–15 Pyroptosis, an inflammasome-mediated pro-inflammatory necrosis, significantly contributes to cartilage degeneration, synovitis, and pain in OA patients.^16,17 Ferroptosis, characterized by iron-dependent lipid peroxidation, and cuproptosis, driven by copper ion-induced aggregation of mitochondria-associated proteins,¹⁸ are discovered novel cell death modalities.¹⁹ Although research on ferroptosis and cuproptosis in OA is limited, epidemiological studies have linked iron and copper overload to an increased risk of OA.²⁰ Notably, inhibiting ferroptosis during the early and middle stages of cartilage injury can restore chondrocyte function and decelerate OA progression. The regulation of cuproptosis plays a significant role in OA diagnosis and subtyping.^20,21 Interestingly, there are interactions between various PCD pathways, such as PANoptosis (a combination of apoptosis, pyroptosis, and necroptosis), which has been implicated in degenerative diseases such as osteoporosis.²² However, there is limited research on the specific relationships between these PCD pathways and OA synovitis, along with the mechanisms underlying their interactions.

Bioinformatics is a multidisciplinary discipline widely used in disease diagnosis, transcriptomics, proteomics, and other fields.²³ Our previous research identified aging-related genes in OA and explored their role in immune inflammation.²⁴ Given the link between diverse PCD pathways and inflammation, autoimmunity reactions in diseases, and their significant role in OA,²⁵ the current study aimed to conduct an in-depth analysis role of the diverse PCD-related genes in OA synovitis and explore their correlation with cluster subgroups, biomarkers, and immune cell infiltration in OA, as outlined in the study’s flowchart (Figure 1). These discoveries could present novel avenues for the early diagnosis and targeted treatment of OA.

Figure 1 Flow Chart. Flowchart for Comprehensive Analysis of PCD Related Genes in Osteoarthritis. *P < 0.05, **P < 0.01, and ***P < 0.001.

Materials and Methods

Data Collection and Processing

We obtained gene expression datasets from the Gene Expression Omnibus (GEO) based on the following criteria: (I) mRNA expression profiles generated through microarray analysis or high-throughput sequencing; (II) tissue samples in the dataset originated from Homo sapiens, including both controls and OA patient groups; and (III) the combined sample size of controls and OA patients in each dataset consisted of at least 10 cases. We obtained seven gene microarray datasets (GSE55235, GSE55457, GSE12021, GSE55584, GSE1919, GSE82107, GSE46750, GSE63359, GSE48556, and GSE32317), one high-throughput RNA sequencing dataset (GSE89408), and two single-cell RNA-seq datasets (GSE152805, GSE176308), as detailed in Table 1. Four synovial tissue datasets (GSE55235, GSE55457, GSE12021, and GSE55584) were merged using the limma package in R to correct the background and normalize each dataset. Batch effects were eliminated using the R package SVA.²⁶ The merged dataset was used for training, and external validation was performed on datasets including OA synovial tissue (GSE89408, GSE1919, GSE82107), synovial cells (GSE46750), peripheral blood leukocytes (GSE63359), peripheral blood mononuclear cells (GSE48556), and early and late OA synovitis (GSE32317). The single-cell RNA dataset GSE176308 from OA synovial tissue contains three samples: GSM5362559 was derived from the painful region of end-stage OA, GSM5362560 was derived from the painful region of early-stage OA, and GSM5362561 was derived from the non-painful region of early-stage OA. We used these samples to represent advanced OA, mid-stage OA, and early OA synovial tissues for the subsequent analyses, respectively.

Table 1 Descriptive Statistics

Download and Collation PCD-Related Genes

PCD-related genes (apoptosis, necroptosis, autophagy, pyroptosis, ferroptosis and cuproptosis) were obtained from FerrDb (http://www.zhounan.org/ferrdb/current/), HADb (http://autophagy.lu/), MSigDB (https://www.gsea-msigdb.org/gsea/msigdb), the KEGG (https://www.genome.jp/kegg/) database. After removing duplicate genes, 1362 genes related to PCD were obtained for analysis. Details are shown in supplementary Table S1.

Identification of Differentially Expressed Genes (DEGs) for OA

We used the limma package of R to analyze differentially expressed genes in Normal and OA synovial samples in the training set, screened for |log₂FC|>1 and FDR<0.05 as the screening criteria.²⁷ The R package ggplot2 visualized volcano maps and gene expression heatmaps. We obtained PCD-related differential genes (PCD-DEGs) by taking intersections of differential genes with PCD-related genes.

Functional Enrichment Analysis

We used the R packages Clusterprofiler and DOSE, with p-value<0.05, as screening criteria for Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Disease Ontology (DO) enrichment analyses. The R package ggplot2 visualized the results. We used the Metascape database to perform cluster enrichment analysis on PCD- DEGs, including GO and KEGG entries, with statistical significance at P-value < 0.05.²⁸

GSVA Enrichment Analysis

GSVA is a non-parametric method for assessing biological pathway enrichment in transcriptomic gene expression datasets.²⁹ We used the “Hallmark.all.v2022.1.Hs.symbols” and “c2.cp. Kegg. symbols” gene sets from the Molecular Signature Database (MSigDB) by downloading them as the reference set. Enriching Normal and OA synovial samples was evaluated by scoring the HALLMARKS and KEGG pathways using the R package GSVA. A significant enrichment was considered to have a p-value < 0.05.

Protein-Protein Interaction Networks (PPI) and Correlation Analysis of PCD-DEGs

Protein interactions of PCD-DEGs were explored using the STRING database (https://string-db.org/), filtered with a confidence score of 0.4, and visualized with Cytoscape (Version 3.9.1) based on Degree scoring.³⁰ The Pearson correlation test was utilized to examine the interaction of PCD-DEGs at the mRNA level. We obtained information on the location of PCD-DEGs genes on chromosomes from the ENSEMBL database and visualized it using the R package RCircos.

Co-Expression Analysis of PCD-DEGs

The GeneMANIA database (http://genemania.org/) is a gene function prediction tool that uses a large set of functional association data to find other genes associated with it based on input genes.³¹ We upload the PCD-DEGs, setting Max resultant genes = 20 and Max resulting attributes = 10.

Identification of Hub PCD-DEGs

We used LASSO regression, SVM-RFE, and Random Forest trees to further screen and identify Hub PCD-DEGs. LASSO regression is helpful for screening variables, handling high-dimensional data, and solving multicollinearity.³² We used the glmnet R package to obtain the optimalλthrough ten cross-validations to identify the PCD-related signature genes. SVM-RFE is an iterative method that uses embedding to identify the genes with the most discriminatory ability.³³ We used the e1071 R package to optimize feature vector models and obtain PCD-related signature genes. The Random Forest algorithm’s recursive Feature Elimination (RFE) is a non-parametric, supervised classification algorithm.³⁴ The R randomForest package was utilized to construct a decision tree that minimizes errors in filtering critical variables. PCD-related signature genes were identified based on genes with relative importance greater than one after tenfold cross-validation. The Venn diagram shows the Hub PCD-DEGs, the intersection of PCD-related signature genes from the three machine-learning methods.

Construction of a Diagnostic Model for OA

We created a nomogram based on Hub PCD-DEGs gene expression using R package rms to improve the clinical applicability of the Hub PCD-DEGs diagnostic prediction OA. We plotted calibration, clinical decision, and clinical impact curves to assess the model’s clinical predictive benefit.

Receiver Operating Characteristic (ROC) Curve Analysis and Differential Expression of Hub PCD-DEGs

We used the R package pROC to perform ROC analyses on the training set and external validation datasets (GSE89408, GSE1919, GSE82107, GSE46750, GSE63359, GSE48556, and GSE32317) to evaluate the diagnostic value of the Hub PCD-DEGs and Model. We used the Wilcox test to analyze the differential expression of Hub PCD-DEGs in synovial samples from Normal and OA patients.

Patients Samples

The human synovial tissues were obtained from 10 OA patients who underwent total knee replacement and relatively normal synovial tissue s from 6 patients with meniscus injury at Guangzhou Red Cross Hospital of Jinan University. All patients signed informed consent, and samples were collected, processed, and analyzed under the guidance of the ethics committee of the Guangzhou Red Cross Hospital of Jinan University (Ethics number 2024–081-01).

qRT-PCR Validation

We extracted total synovial tissue RNA using Trizol (Servicebio) and then reverse-transcribed the total RNA into complementary DNA (cDNA) using Takara Prime Script ^®RT Master Mix. qRT-PCR was performed using 2× SYBR Green qPCR Master Mix (No ROX) (Service) and using the Easy Cycler PCR instrument of ANALYTIKJENA, according to the manufacturer’s instructions. The GAPDH gene was used as an internal reference gene. Each biological sample was technically repeated three times. The primer sequences are as follows:

GAPDH Forward:5′-ACACCCACTCCTCCACCTTT-3′; Reverse:5′-TTACTACTTGGAGGCCATGT-3′, TNFAIP3; Forward:5′- ACCAGGACTTGGGACTTTGC-3′; Reverse:5′-TGTGCTCTCCAACACCTCTC-3′, PPP1R15A; Forward:5′-TGAGACTCCCCTAAAGGCCA-3′; Reverse:5′-CCAGACAGCCAGGAAATGGA-3′, INHBB Forward:5′-GCGCGTTTCCGAAATCATCA-3′; Reverse:5′-TTCTGGTTGCCTTCGTTGGA-3′, DDIT4 Forward:5′-GGTTCGCACACCCATTCAAG-3′; Reverse:5′-CAGGGCGTTTGCTGATGAAC-3′, JUN Forward:5′-AACAGGTGGCACAGCTTAAAC-3′; Reverse:5′-CAACTGCTGCGTTAGCATGAG-3′.

Construction of a miRNA-TF-mRNA Regulatory Network for Hub PCD-DEGs

We used the miRTarBase,³⁵ Starbase,³⁶ and Targetscan databases³⁷ to predict miRNAs with regulatory relationships to Hub PCD-DEGs, taking the intersection of the three databases. Prediction of transcription factors (TFs) for possible roles of Hub PCD-DEGs using the Enrichr database (https://maayanlab.cloud/Enrichr/),³⁸ with a p-value ≤ 0.05 as a screening condition. The miRNA-TF-mRNA regulatory network was visualized using Cytoscape (3.9.1) sorted by Degree rank.

GSEA Enrichment Analysis

Gene Set Enrichment Analysis (GSEA) is a computational method to identify statistically significant differences in gene sets between two biological states.³⁹ We divided the training set into a high-expression group (≥50%) and a low-expression group (<50%) based on the expression of Hub PCD-DEGs. The “Hallmark.all.v2022.1.Hs.symbols” reference set from the MSigDB database was used to evaluate pathways and molecular mechanisms associated with differential Hub PCD-DEGs expression between groups.

Identification of PCD Clusters in OA

We performed unsupervised clustering of OA synovial samples using the ConsensusClusterPlus R package based on the expression values of PCD-DEGs with K set to 10. The optimal K value is determined by analyzing the full cumulative distribution function (CDF) index. The distribution of samples between clusters is evaluated through principal component analysis (PCA) and the expression of PCD-DEGs in two PCD Clusters using the Wilcox test.

PCD Score Construction

We created a scoring system using PCA to describe PCD Clusters in OA patients based on PCD-DEGs.

PC1 indicates principal component 1, PC2 indicates central component 2, and I reflect the expression of subgroup-specific genes.

Molecular Characterisation and Weighted Gene Co-Expression Network Analysis (WGCNA) in PCD Clusters

We used reference gene sets from MSigDB, including “c5.go.v7.5.1.entrez.gmt”, “c2.cp.Kegg.symbols.gmt”, and “Hallmark.all.v2022.1.Hs.symbols.gmt”, for GSEA enrichment analysis to explore two PCD clusters and associated biological functions. We constructed a weighted co-expression network using the R package WGCNA and analyzed genes with an absolute median difference in expression greater than 0.25 across different PCD Clusters.⁴⁰ The genes were grouped into other modules according to their TOM values. The module most relevant to the cluster was chosen, and the genes in that module were screened for Hub genes using cor.MM>0.7 and cor.GS>0.5 as criteria. The module’s molecular biology of the Hub genes was analyzed using GO and KEGG enrichment.

Immune Characteristics of OA and PCD Clusters

Single sample gene set enrichment analysis (ssGSEA) is widely used in immuno-infiltration-related bioinformatics studies.⁴¹ We used the R package GSVA to calculate enrichment scores for immune cells and functions in Normal and OA synovial samples and PCD Clusters. The relative abundance of immunocytes, the enrichment scores of immune functions, the status of the HLA gene, and chemokine in Normal and OA samples were examined by the Wilcoxon test. ESTIMATE assessed three scores: Stromal, Immune, and Immune Microenvironment Score for patients with OA and PCD Clusters.⁴² The Spearman test was applied to analyze the correlation between Hub PCD-DEGs and immune cells and immune function.

Single-Cell RNA Sequencing Analysis

We used the Seurat R package to analyze single-cell RNA datasets of OA synovial tissues (GSE152805 and GSE176308). Data underwent quality control using the following screening criteria: (I) cells with gene expression in ≤3 cells; (II) cells expressing <200 genes or >2500 genes; and (III) cells with ≥5% mitochondria and ribosomes. The data were normalized using the NormalizeData function and highly variable genes in the cells were identified based on the mean-variance relationship. The cells were clustered based on the first 2000 variable genes using PCA (resolution=0.6) pairs. The genes underwent differential analysis using FindAllMarkers, screening with min. pct = 0.25 and logfc.threshold = 0.25. Diverse cell clusters were manually identified and annotated using the CellMarker database (http://biocc.hrbmu.edu.cn/CellMarker/), combined with relevant literature, based on each cluster’s top ten most up-regulated genes. The proportion of cell clusters and the expression of Hub PCD-DEGs in different samples and cell populations were also calculated and visualized.

Statistical Analysis

All statistical analyses were performed using R version 4.2.2. The Wilcoxon test was utilized to compare the two groups statistically. The qRT-PCR statistical analyses were presented as the mean ± standard deviation for a minimum of three individual experiments. The statistical significance of differences was determined through the unpaired, two-tailed Student’s t-test. Statistical significance was indicated by *P < 0.05, **P < 0.01, and ***P < 0.001.

Results

Data Processing and Differential Gene Screening

We combined four synovial datasets from platform GPL96 as the training set, which included 65 samples: 29 Normal and 36 OA synovial samples. UMAP analysis showed that the batch effect of the merged dataset was well eliminated (Figure 2A, and B). We used the R package limma to identify 278 differentially expressed genes, with 108 up-regulated and 170 down-regulated in OA. The results were visualized by volcano plots (Figure 2C), where the top 30 genes with the most significant differences were presented in the heat map (Figure 2D).

Figure 2 Data processing and DEGs screening. (A and B) Diagram of the distribution of each sample before and after removal of the batch effect. (C) Heat map of the top 30 DEGs, with the left half representing Normal synovial samples, the right half representing OA synovial samples, and red and blue representing upregulation and downregulation, respectively. (D) Volcano map of DEGs. Red nodes indicate up-regulated DEGs, blue nodes indicate down-regulated DEGs and grey nodes indicate genes that are not significantly differentially expressed.

Enrichment Analysis of OA

We performed GO, KEGG, and DO enrichment analyses on DEGs. GO enrichment analysis showed that the differential genes were mainly enriched in fat cell differentiation, collagen-containing extracellular matrix, RNA polymerase II transcription regulator complex, cytokine activity, signaling receptor activator activity, and significant receptor activity, cytokine activity, and signaling receptor activator activity (Figure 3A). KEGG enrichment analysis showed that DEGs were mainly enriched in various immune-inflammatory pathways, such as the TNF signaling pathway, IL-17 signaling pathway, Rheumatoid arthritis, Osteoclast differentiation, AGE-RAGE signaling pathway in diabetic complications, NF-kappa B signaling pathway, and Cytokine-cytokine receptor interaction (Figure 3B and C). DO enrichment results suggested that DGEs were associated with diseases such as myeloid neoplasm, ulcerative colitis, bone marrow cancer, osteoarthritis, and colitis (Figure 3D). We used GSVA analysis to explore the enrichment of HALLMARKS and KEGG pathways in OA. Compared with the Normal group, the enrichment of the HALLMARKS pathway showed that TNFA signaling via the NFKB pathway, apoptosis, and P53 pathway was significantly upregulated in OA. In contrast, DNA repair, oxidative phosphorylation, and bile acid metabolism were downregulated considerably (Figure 3E). KEGG pathway scoring showed that small cell lung cancer, adipocytokine signaling pathway, and chronic myeloid leukaemia were significantly upregulated in OA. In contrast, several metabolic pathways were downregulated in OA, such as butanoate metabolism, starch and sucrose metabolism, selenoamino acid metabolism, starch and sucrose metabolism, and selenoamino acid metabolism (Figure 3F).

Figure 3 Enrichment Analysis of OA. (A) The size of the bubbles represents the number of DEGs (the more significant the circle, the more DEGs), the color represents the p-value (the redder the shade, the more significant the p-value described), and the circle represents biological processes (BP), the square represents cellular components (CC). The triangle represents molecular function (MF). (B) KEGG enrichment analysis bubble map showing the top 10 most considerable pathway enrichments. (C) KEGG pathway interactions map. (D) DO enrichment analysis bubble plots. (E) Bar map of GSVA enrichment scoring of HALLMARKS pathways in OA patients relative to Normal human synovial tissue: red represents pathways upregulated in OA, and blue represents downregulation. The higher the enrichment scoring, the longer the bar. (F) Scored bar graph of KEGG pathway GSVA enrichment.

PPI and Correlation Analysis of PCD-DEGs

We took the intersection of DEGs with PCD-related genes and obtained 26 PCD-DEGs, of which four were up-regulated in OA and 22 were down-regulated in OA (Figure 4A). The chromosomal location of PCD-DEGs is shown in Figure 4B. We conducted a PPI analysis on PCD-DEGs, ranked by Degree scores and visualized using Cytoscape. The analysis revealed that IL6, IL1B, and JUN had the most robust interactions with other proteins in the network (Figure 4C). We analyzed PCD-DEGs at the mRNA interaction level using the Spearman correlation test. The results showed that PCD-DEGs acted tightly at the protein and mRNA levels, with MCL1 significantly positively correlated with ZFP36 (r = 0.88, P = 5.04e-12) and DUOX2 significantly negatively correlated with CX3CR1 (r = −0.66, P =2.01e-05) (Figure 4D).

Figure 4 PPI and Correlation Analysis of PCD-DEGs. (A) Heat map of differential expression of PCD-DEGs, * indicates p < 0.05, ** indicates p < 0.01; *** represents p < 0.001. (B) Chromosome location map of PCD-DEGs. (C) PPI analysis of PCD-DEGs, the larger the gene node, the darker the red, the higher the gene Degree score, and the stronger the reciprocal relationship. (D) Correlation study of PCD-DEGs.

Enrichment Analysis of PCD-DEGs

We explored the molecular mechanism of PCD-DEGs in OA by GO and KEGG enrichment analysis. KEGG enrichment analysis showed that the enrichment of PCD-DEGs was consistent with the enrichment of DEGs for rheumatoid arthritis, IL-17 signaling pathway, TNF signaling pathway, NF-kappa B signaling pathway, and osteoclast differentiation (Figure 5A). GO enrichment results showed that PCD-DEGs are mainly involved in the regulation of the apoptotic signaling pathway, cellular response to external stimulus, neuroinflammatory response, response to hypoxia, intrinsic apoptotic signaling pathway, cytokine activity, signaling receptor activator activity, cytokine receptor binding and other biological functions (Figure 5B and C). To further reveal the biological properties of PCD-DEGs, based on GeneMANIA, we constructed and analyzed the network of PCD-DEGs and their co-expressed genes. PCD-DEGs exhibit a complex web of PPI interactions: co-expression, 83.57%; physical interaction, 10.34%; co-localization, 4.67%; genetic interaction, 1.13%; Prediction, 0.26%; pathway, 0.01% (Figure 5D).

Figure 5 Enrichment Analysis of PCD-DEGs. (A) KEGG chord diagram of PCD-DEGs. (B) GO chord diagram of PCD-DEGs. (C) GO, KEGG clustering analysis of PCD-DEGs. (D) PCD-DEGs and their co-expression genes were analyzed via GeneMANIA.

Identification of Hub PCD-DEGs

Three machine-learning algorithms were applied to screen for hub PCD-DEGs. The ten most minor binomial deviation PCD signature genes were identified by lasso regression (Figure 6A and B). The SVM-RFE results showed the highest accuracy (0.971) and the lowest error (0.0286) in model effectiveness when selecting the top 24 (Figure 6C and D). The RF was used to rank the importance of the PCD-DEGs (Figure 6E and F), and 12 PCD signature genes were obtained with a volume greater than one as the screening condition. Finally, five intersected Hub PCD-DEGs, namely TNFAIP3, JUN, PPP1R15A, INHBB, and DDIT4, were obtained through three algorithms (Figure 6G).

Figure 6 Machine-learning Algorithms to Screen for Hub PCD-DEGs. (A and B) Construction of PCD signature genes using LASSO regression. (C and D) The plot of maximum accuracy and minimum error of the SVM-RFE algorithm for screening optimal parameters. (E) The relationship between the number of random forest trees and the error rate. (F) Relative importance scores for PCD-DEGs. (G) LASSO, Random Forest, and SVM-RFE algorithms for screening Venn diagrams of PCD feature gene intersections.

Construction of a Diagnostic Model and ROC Curve for OA

We constructed a diagnostic Model for OA based on the expression of five Hub PCD-DEGs (Figure 7A). The clinical calibration curve showed a high level of agreement between the actual and predicted risks of the Model (Figure 7B). The results of the clinical decision curves revealed that the Model has a high net benefit as the threshold probability increases (Figure 7C). The Clinical Impact Curve (CICA) results showed that the model predicts a minimal difference between the number of people with OA and the actual number of people with OA (Figure 7D). These results indicated that the Model constructed using Hub PCD-DEGs has high clinical predictive efficacy for OA. We used ROC curve analysis to evaluate further the diagnostic performance of Hub PCD-DEGs and Model in the training and three validation sets. The results indicate high diagnostic value for TNFAIP3, JUN, PPP1R15A, INHBB, DDIT4, and nomogram in both the training and validation sets, with AUC more significant than 0.6 for all five Hub PCD-DEGs and nomogram (Figure 7E-G and Supplementary Figure S1A). In order to delve deeper into the capacity of Hub PCD-DEGs to distinguish between early and late-stage OA, we developed a predictive model using the GSE32317 dataset (Figure 7H), and corroborated the model’s accuracy with calibration plots (Figure 7I) and ROC curves (Figure 7J).The results indicated that TNFAIP3, JUN, PPP1R15A, INHBB, and Model had high discriminatory ability. The model had the highest value, with an AUC of 0.95. In addition, we performed external validation on the data sets GSE46750, GSE48556, and GSE63359, and the results showed that Hub PCD-DEGs and diagnostic models have predictive value in synoviocytes cultured in vitro and peripheral blood of OA patients (Supplementary Figure S1B-D).

Figure 7 Construction of a Diagnostic Model and ROC Curve for OA. (A) Nomogram for diagnosing OA. (B) Calibration curve to assess the predictive accuracy of the Nomogram (the closer to the ideal dashed line, the more reliable the result). (C) The clinical decision curve of the diagnostic model (the further the red line endpoints are from the grey line, the higher the net benefit of the model). (D) The clinical impact curve of the diagnostic model (the more significant the area under the curve, the more reliable the diagnostic model). (E) ROC curve analysis of Hub PCD-DEGs and Nomogram in training sets. (F and G) ROC and clinical decision curve analysis of Hub PCD-DEGs and Nomogram in GSE89408 and GSE1919 datasets. (H) Nomogram for identifying early and late stages of OA in GSE32317 datasets. (I) Identifying early and late calibration curves for OA using a nomogram. (J) ROC curve analysis of Hub PCD-DEGs and Nomogram in GSE32317 datasets.

Differential Expression Analysis and qRT-PCR Validation of Hub PCD-DEGs

We further analyzed the differences in the expression of Hub PCD-DEGs between Normal and synovial samples in the training and validation sets using the Wilcoxon test, which showed that Hub PCD-DEGs were significantly downregulated in the training set and two validation sets of OA synovial samples (Figure 8A–C). PPP1R15A and INHBB were downregulated in late OA synovial samples compared to early samples (Figure 8D). We performed qRT-PCR on total mRNA from synovial membranes of six patients with meniscal injuries and ten patients with OA to further validate the mRNA expression levels of Hub PCD-DEGs. The results showed that all Hub PCD-DEGs were significantly downregulated in the OA synovial samples, consistent with the expression in the training set and the two validation sets (Figure 8E).

Figure 8 Differential Expression Analysis and qRT-PCR Validation of Hub PCD-DEGs. (A-C) Differential expression of the Hub PCD-DEGs in the training set, GSE1919, and GSE89408 datasets. (D) Differential expression of the Hub PCD-DEGs in the early and end-stage OA. (E) The qRT-PCR method detected the mRNA expression levels of five Hub PCD-DEGs. *P < 0.05, **P < 0.01, and ***P < 0.001.

GSEA Enrichment Analysis of Hub PCD-DEGs

To better understand the potential biological functions of Hub PCD-DEGs, we conducted GSEA analysis separately. Our results showed that the downregulated groups of TNFAIP3, JUN, PPP1R15A, INHBB, and DDIT4 were significantly enriched in different immune responses and inflammatory pathways, including Inflammatory Response, Interferon Alpha Response, Complement, and others (Figure 9A–E).

Figure 9 GSEA Enrichment Analysis of Hub PCD-DEGs. (A-E) Ridge maps of TNFAIP3, JUN, PPP1R15A, INHBB, DDIT4 GSEA HALLMARK pathway enrichment, respectively.

Predicting miRNA and TF of Hub PCD-DEGs

We predicted the miRNAs and TFs that Hub PCD-DEGs might act on using an online database to visualize the regulatory network using Cytoscape. The network contains 82 nodes involving 204 edges, 64 miRNAs, and 14 TFs. hsa-miR-26a-5p, hsa-miR-1297, hsa-miR-149-5p, hsa-miR-495-3p, and hsa-miR-26b-5p may be the acting target miRNAs co-regulated by the Hub PCD-DEGs (Figure 10).

Figure 10 miRNA-TF-mRNA Regulatory Network of Hub PCD-DEGs. Red circles represent genes; blue V-shapes represent predicted miRNAs; green diamonds represent TFs.

Immune Characteristics of OA

We used ssGSEA analysis to assess the level of infiltration of 28 immune cell species and activation of 13 immune functions in the training set and the GSE89408 dataset. The results showed that compared to Normal synovial samples, CD56bright natural killer cell, Gamma delta T cell, Immature dendritic cell, Macrophage, Natural killer cell, Regulatory T cell, and Type 1 T helper cell were upregulated in OA. In contrast, Central memory CD8 T cell was downregulated in OA. HLA and T cell co-inhibition were significantly activated in OA (Figure 11A and B). We utilized Spearman correlation to analyze the relationship between immune cells and immune function. Our findings indicated significant correlations between various immune cells and immune function. For instance, there was a notable positive correlation between macrophages and MDSC (r = 0.87, p= 6.33e-12). Additionally, HLA was positively associated with Immature B cell (r = 0.80, p= 6.77e-09), please refer to Supplementary Figure S2A, C. We further assessed OA immune scores using ESTIMATE analysis, which showed that OA synovial samples had higher Stromal, immune, and microenvironmental scores than Normal, please refer to Supplementary Figure S2B, D. We used Spearman correlation analysis to investigate the correlation between Hub PCD-DEGs and immune cells and immune function. The results showed that Hub PCD-DEGs were significantly and negatively correlated with various immune cells and immune functions, such as macrophages (Figure 11C and D). Higher expression levels of HLA genes, such as HLA-DQB1, HLA-DMA, and HLA-DRA (Figure 11E and F), and higher expression levels of chemokine, such as CXCL6 (Figure 11G and H), were detected in OA compared to Normal controls.

Figure 11 ssGSEA Immune Infiltration Analysis. (A and B) Violin plots of differences in the infiltration of 28 immune cells and 13 immune functions in Training Sets and the GSE89408 dataset. (C and D) Correlation analysis of Hub PCD-DEGs with 28 immune cells and 13 immune functions in Training Sets and the GSE89408 dataset, * denotes p < 0.05, ** denotes p < 0.01, *** indicates p < 0.001. (E and F) The expression differences of each HLA gene in Training Sets and the GSE89408 dataset. (G and H) The expression differences of each chemokine gene in Training Sets and the GSE89408 dataset.

Identification and Molecular Characterization of PCD Clusters

We used unsupervised cluster analysis to divide the OA synovial samples into two PCD Clusters: A with high expression (n=16) and B with low expression (n=20) (Figure 12A and B). PCA analysis showed that the two Clusters had distinctly different expression patterns (Figure 12C). We further visualized the expression levels of PCD-DEGs in both Clusters using box plots showing that Cluster A had higher PCD-related gene expression than Cluster B (Figure 12D). However, Cluster B has a higher PCD score than Cluster A (Figure 12E). To compare the molecular characteristics of two PCD Clusters, we utilized GSEA enrichment to evaluate HALLMARKS, KEGG pathway, and GO enrichment between them. Hallmark pathway enrichment results showed that allograft rejection, interferon alpha response, interferon-gamma response, and inflammatory response were upregulated in Cluster B (Figure 12F). KEGG pathway enrichment results indicate that the intestinal immune network for IgA production primary immunodeficiency, allograft rejection, and asthma are upregulated in Cluster B (Figure 12G). GO functional enrichment results indicated that the B cell receptor signaling pathway, B cell-mediated immunity, Adaptive immune response, leukocyte-mediated immunity, and Immune response regulating signaling pathway were upregulated in Cluster B (Figure 12H).

Figure 12 Identification of PCD Clusters. (A) Relative changes in the consensus clustered cumulative distribution function (CDF) and area under the curve for k = 2–10. (B) Heat map of an OA synovial sample with k = 2. (C) PCA analysis of two Clusters. (D) Box plots of the differential expression of PCD-DEGs between the two Clusters. *P < 0.05, **P < 0.01, and ***P < 0.001. (E) Differences in PCD scoring for the two Clusters. GSEA enrichment analysis of PCD Clusters. (F) Ridge plot of GSEA enrichment differences in the HALLMARKS pathway for PCD Clusters, with enrichment scores >0 indicating upregulation in Cluster B, and <0 downregulation in Cluster B. (G) Differential ridge map of GSEA enrichment in the KEGG pathway of the PCD Clusters. (H) Top 5 GO entries enriched by GSEA in Cluster A and Cluster B, respectively.

WGCNA Analysis of PCD Clusters

To further explore the role of PCD Clusters in OA, we used the R package WGCNA to identify the most relevant gene modules for Clusters. The results showed that the turquoise module (1437 genes) had the strongest correlation with cluster (cor = 0.55; P= 6e-04) (Figure 13A–D). We selected genes with cor.MM>0.7 and cor.GS>0.5 in the turquoise module as the Hub genes in the Key module (Figure 13E). KEGG enrichment analysis showed that the Hub genes were significantly enriched in Osteoclast differentiation, B cell receptor signaling pathway, Rheumatoid arthritis, Phagosome, and Fc epsilon RI signaling pathway (Figure 13F). GO enrichment analysis showed that Hub genes were significantly enriched in BPs such as myeloid cell activation involved in the immune response, myeloid leukocyte activation, leukocyte activation involved in the immune response, and cell activation involved in the immune response. Cellular Components (CC) include tertiary granules, endocytic vesicles, and secretory granule membranes. In Molecular Functions (MF), eg actin binding, inhibitory MHC class I receptor activity, and MHC class I receptor activity (Figure 13G).

Figure 13 WGCNA analysis of PCD Clusters. (A) Sample dendrogram and trait heatmap. (B) Determine the optimal soft threshold and construct a scale-free network by choosing nine based on where the red line (R²0.9) is located. (C) MAD > 0.25 in the top 25% of the gene cluster dendrogram. (D) Heat map of module-trait relationships, where each color represents a co-expression module, where values indicate module-trait correlation coefficients and p-values. (E) Scatter plot of correlations between gene significance (GS) and module membership (MM) in the turquoise module. (F) Chordal map of KEGG enrichment of crucial genes within the turquoise module. (G) Histogram of GO enrichment of critical genes within the turquoise module.

Immune Characteristics of PCD Clusters

We used ssGSEA analysis to assess the differences in immune cell infiltration and immune function activation between the two PCD Clusters. The results showed that Activated B cells, Activated CD8 T cells, Immature B cells, Immature dendritic cells, MDSC, Macrophage, Regulatory T cells, Effector memory CD4 T cell, and Memory B cell significantly infiltrated in Cluster B. Check-point, HLA, Inflammation-promoting, Parainflammation, T cell co-stimulation immune functions were higher at Cluster B activation levels (Figure 14A). The ESTIMATE score suggests that Cluster B has a higher immunological and microenvironmental score than Cluster A (Figure 14B). Higher expression levels of HLA genes, such as HLA-DQA1, HLA-DMA, and HLA-DOA (Figure 14C) and higher expression levels of chemokine, such as CXCL5, CXCL10 and CXCL22 (Figure 14D and E), were detected in Cluster B compared to Cluster A.

Figure 14 Immune infiltration analysis of ssGSEA for PCD Clusters. (A) Heat map of differences in the distribution of immune cells and immune function in PCD Clusters. (B) ESTIMATE Scores for the PCD Clusters. (C) Violin plots of differences in immune cells and immune function between different PCD Clusters. (D) The expression differences of each HLA gene in PCD Clusters. (E) The expression differences of each chemokine gene in PCD Clusters.

Single-Cell RNA Sequencing Analysis of the OA Synovium

The expression profiles of 3688 genes were obtained after quality control of the GSE152805 single-cell RNA dataset (Figure 15A). For detailed quality control results of the GSE152805 dataset, please refer to Supplementary Figures S3A-C. PCA identified nine cell subpopulations: Monocytes, Epithelial cells, Chondrocytes, Macrophages, Natural killer T cells, Tissue stem cells, CD4+ T cells, B cells, and NK cells (Figure 15B and C). Notably, we found differential expression of Hub PCD-DEGs in different cell clusters of the OA synovium (Figure 15D). JUN, PPP1R15A, and DDIT4 genes were highly expressed in different cell clusters. In contrast, INHBB was barely expressed in all cell types in OA. TNFAIP3 was predominantly expressed in macrophages and B cells. To better understand the distribution of cell clusters and the differences in the expression of Hub PCD-DEGs during various stages of OA, we performed single-cell RNA analysis on the GSE176308 dataset. For detailed quality control results of the GSE176308 dataset, please refer to Supplementary Figure S4A-D. The synovium in advanced, mid-stage and early OA contained macrophages, monocytes, DCs, Stromal cells, NKT cells, and mast cells. The synovial cell composition varies dramatically at different stages of OA (Figure 15E and G). As shown in Figure 15F, macrophages dominate early OA synovium, and the proportion of various cells tends to even out with the progression of OA. In contrast, advanced OA synovium exhibits a predominance of Stromal cells. Hub PCD-DEGs were inconsistently expressed in cell clusters at different stages of OA. As shown in Figure 15H and Supplementary Figure S3E, the expression of TNFAIP3 and DDIT4 was significantly decreased in advanced synovial samples.

Figure 15 Single-cell RNA sequencing analysis of the OA synovium. (A) Heatmap of the top 10 marker gene expression in different cell clusters of the GSE152805 dataset. (B and C) Cells in the GSE152805 dataset were clustered into nine types via the t-SNE dimensionality reduction algorithm. Each color represented the annotated phenotype of each cluster. (D) Feature and violin plots showing the distribution of Hub PCD-DEGs in various cell types in the GSE152805 dataset. (E) Cells in the GSE176308 dataset were clustered into six types via the t-SNE dimensionality reduction algorithm. (F and G) Proportion and distribution of cell clusters in OA synovial samples from different stages. (H) Feature plots showing the distribution of Hub PCD-DEGs in various cell types in the GSE176308 dataset.

Discussion

OA is a common cause of chronic pain and disability in older people. The complexity of its pathogenesis and the absence of reliable methods for diagnosing, staging, and monitoring joint pathological changes remain significant barriers to advancing OA treatments.⁴² Emerging research suggests that diverse PCD pathways play a critical role in OA progression.⁴³ However, comprehensive investigations into the diverse PCD pathways in OA are still limited.

In this study, we analyzed the Normal and OA synovial gene expression datasets from the GEO database and identified 278 DEGs and 26 PCD-DEGs. GO and KEGG enrichment analyses revealed that DEGs and PCD-DEGs were enriched for immune-inflammatory pathways, including the TNF signaling pathway, IL-17 signaling pathway, NF-kappa B signaling pathway, and cytokine activity, which is consistent with previous studies in OA.⁴⁴ TNF-α plays a key role in regulating inflammatory immune responses, with activation of the TNF signaling pathway driving OA-related inflammation through diverse PCD pathways.⁴⁵ The activation of the IL-17 signaling pathway stimulates chemokine release from chondrocytes and synovial fibroblasts, contributing to cartilage degradation and immune cell infiltration in the synovium.⁴⁶ Furthermore, IL-17 signaling induces mitochondrial dysfunction in fibroblast-like synoviocytes (FLS), leading to autophagic apoptosis.⁴⁷ NF-κB signaling has been implicated in synovial inflammation and chondrocyte apoptosis during OA progression.⁴⁸ Blocking this pathway can reduce chondrocyte apoptosis and pyroptosis by inhibiting the activation of NLRP3 inflammasomes.⁴⁹ These findings suggest that the transcription of PCD-DEGs may play a role in the development of synovial inflammation in OA.

We used three machine learning algorithms and ROC curve analyses to screen and validate TNFAIP3, JUN, PPP1R15A, INHBB, and DDIT4 as novel biomarkers for the diagnosis of OA. Tumor Necrosis Factor Alpha-Induced Protein 3 (TNFAIP3), a zinc finger protein and ubiquitin editing enzyme, inhibits NF-κB activation and apoptosis induced by cytokines such as TNF-α and IL-1β or via Toll-like receptors (TLRs) and is involved in immune and inflammatory responses.⁵⁰ It was found that TNFAIP3 knockout mice exhibited basal and LPS-induced higher levels of Nlrp3 expression, which exacerbated synovial macrophage pyroptosis and RA-associated inflammatory response and cartilage destruction through activation of the Nlrp3 inflammasome/caspase-1/IL-1 signaling axis.⁵¹ Jun Proto-Oncogene (JUN) is a member of the nuclear transcription-activating protein l family, which is involved in various physiological processes such as cell cycle progression, differentiation, and apoptosis.⁵² JUN and BATF bind to form a complex that regulates anabolic and catabolic gene expression in chondrocytes, contributing to OA cartilage destruction.⁵³ JUN overexpression was also associated with the inhibition of SOX9 transcriptional activity and reduced expression of type II collagen in chondrocytes.⁵⁴ Protein phosphatase one regulatory subunit 15A (PPP1R15A) is a member of a protein family that modulates the expression of stress-responsive genes by growth arrest and DNA damage signals.⁵⁵ PPP1R15A plays a crucial role in lysosomal biogenesis and maintaining autophagic flux during starvation.⁵⁶ Additionally, activation of PPP1R15A enhances autophagy and inhibits apoptosis stimulated by LPS combined with amino acid deprivation by regulating the mTOR signaling pathway in macrophages.⁵⁷ In RA, overexpression of PPP1R15A is associated with circulating anti-citrullinated protein antibodies (p = 0.030), which promote inflammatory responses.⁵⁸ Inhibin subunit beta B (INHBB), a member of the transforming growth factor beta (TGF-β) superfamily, is involved in the synthesis of inhibitors and activators, mainly in endocrine, reproductive, apoptosis and tumor processes, and is considered to be a valuable biomarker for various cancer types.⁵⁹ Overexpression of KLF10 may inhibit chondrocyte proliferation and migration in mice by suppressing INHBB expression and is involved in the progression of OA.⁶⁰ DNA Damage Inducible Transcript 4 (DDIT4) is a stress regulator protein induced by hypoxia, endoplasmic reticulum stress, cellular iron depletion, and DNA damage.⁶¹ DDIT4 is also thought to be an endogenous regulator of mTOR signaling and autophagy in articular chondrocytes.⁶² Senescent and OA articular cartilage, synovium, and meniscus exhibit downregulation of DDIT4. DDIT4 knockout mice show higher synovitis scores, cartilage destruction, and bone redundancy production.⁶³ Interestingly, all five diagnostic biomarkers were significantly down-regulated in OA synovial samples by the Wilcox test, which was also verified by experimental qRT-PCR. GSEA analysis to find that the Hub PCD-DEGs down-regulated group exhibited enrichment in multiple immune responses and inflammatory pathways. However, whether their down-regulation is involved in the development of OA synovial inflammation needs further experimental studies. Our study investigated the molecular mechanisms underlying PCD-related diagnostic in OA and their regulatory network with miRNA and transcription factors. This information can guide the development of targeted therapies for OA and immunotherapy.

Although OA is not an immune-related disease, immunity plays a pathogenic role in its development. Several studies have confirmed that the lining beneath the synovium in patients with OA is rich in macrophages, T cells, and small numbers of mast cells, B cells, and plasma cells.^6,64 Our findings also support this, as we observed higher levels of infiltration of macrophages, Activated CD8 T cells, Regulatory T cells, Activated B cells, and Memory B cells in OA synovial samples. Macrophages, the predominant innate immune cells in the synovium of OA, are present in the knee joints of 76% of OA patients and significantly correlate with knee pain and OA imaging severity and osteoclastogenesis.^65,66 CD4 + T cells dominate T cell infiltration in synovial tissue, and activation of CD4 + T cells induces the expression of macrophage inflammatory protein-1γ (MIP-1γ) and NF-κB, leading to increased macrophage infiltration and osteoclast formation, which exacerbates articular cartilage damage.⁶⁷ CD8 + T cells are thought to induce the expression of tissue inhibitor of metalloprotein-1 (TIMP-1) in correlation with the severity of OA.⁶⁸ Analysis of B cells in OA synovium shows an oligoclonal nature, indicating antigen-driven activation.⁶⁹ This suggests that B cells may contribute to OA development and progression rather than being bystanders. We found that 5 hub PCD-DEGs had a significant negative correlation with immune cells and functions, further implying potential as a new target to prevent OA immune inflammation.

To investigate the mechanism of the role of PCD in OA, we identified two PCD-related Clusters of OA by unsupervised consensus clustering analysis. GSEA enrichment analysis revealed that Cluster B exhibited multiple adaptive immune responses and inflammation-related pathways upregulated in Cluster B, such as B cell receptor signaling pathway, B cell-mediated immunity, Adaptive immune response, Leukocyte mediated immunity, Immune response regulating signaling pathway, and inflammatory response. The hub genes in the co-expressed gene module highly associated with Cluster B were significantly enriched in multiple immune reactions, consistent with GSEA enrichment analyses, which suggest that adaptive immune responses may play an essential role in the pathogenic mechanisms of Cluster B. Several studies have found that serum autoantibodies against chondrocyte surface proteins have been detected in the blood of OA patients using ELISA, further confirming the activation of humoral adaptive immune responses in OA.^70,71 We further analyzed the immune microenvironmental characteristics of PCD Clusters using the ssGSEA algorithm, which also confirmed the above results. The molecular differences between the two PCD Clusters require further investigation through experimental studies.

Single-cell RNA sequencing analysis has unveiled a rich diversity of cell types within the synovium of individuals with OA, including Monocytes, Epithelial cells, Macrophages, Stromal cells, and Mast cells. This discovery reinforces the findings from prior research and aligns with the results obtained from our ssGSEA immune cell infiltration analysis.⁷² Interestingly, we found that macrophages were mainly infiltrated in the synovial membrane in the early stages of OA. However, as OA progresses, various other cell types in the advanced OA synovium undergo a gradual replacement by synovial stromal cells. This indicates the pathological process in which OA synovial histology gradually changes from cell proliferation and immune cell infiltration in the synovial lining to interstitial vascularization and fibrosis under the synovial lining.⁷³ This transformation may be related to the disease’s duration and metabolic and structural changes in other joint tissues.⁷⁴ Notably, Hub PCD-DEGs displayed varying expression patterns in different synoviocyte populations and at different stages of OA, which deviated from the results obtained through transcriptome sequencing. Therefore, it is essential to consider the variability among different cell clusters to comprehensively understand the mechanisms underlying OA development and progression when investigating OA’s pathogenesis and the role of Hub PCD-DEGs.

However, there are limitations to this study. The transcriptomic data from this study were obtained from public databases that do not allow for more clinically relevant information, and the potential influence of heterogeneity in patient populations and clinical characteristics on the results of this study cannot be excluded. The small sample size analyzed compromised the accuracy of the results, and a study with a large sample size and prospective nature needs to be designed to validate the model. We used synovial tissue from patients with meniscal tears as controls to represent the “physiological state of normal synovial tissue”. However, it is important to note that meniscal tears are a known risk factor for OA, and synovial inflammation is reported to be present in up to 50% of patients with meniscal tears.^75,76 This could introduce potential bias into the experimental results. We only validated through database analysis and simple experiments, lacking support from further in vivo and in vitro experiments.

Conclusion

Our study revealed a connection between PCD-related genes and immune inflammation in OA synovial tissue. We identified Hub PCD-DEGs, including TNFAIP3, JUN, PPP1R15A, INHBB, and DDIT4, which may serve as emerging early diagnosis and treatment targets for OA. Additionally, we characterized PCD clusters associated with the OA immunoinflammatory phenotype and uncovered dynamic cellular changes within OA synovial tissue, opening up the possibility of precisely targeted therapy for OA patients. However, further studies are needed to support our findings.