Microarray datasets for RA were retrieved from the GEO database using the following criteria: “Rheumatoid Arthritis,” “Expression profiling by array,” “Homo sapiens,” and “sample count >20.” Three gene expression profiles (GSE12021, GSE55457, and GSE55235) based on the GPL96 platform were selected for further analysis. GSE12021 includes 12 synovial tissue samples from RA patients and 9 from healthy knee joints. GSE55457 consists of 23 samples, with 10 from healthy synovial tissue and 13 from RA patients. GSE55235 contains 20 samples, 10 from healthy knee joints and 10 from RA patients. Furthermore, GSE77298 was used as an independent dataset because the validation group (including RA/normal synovial samples) consisted of 16/7 samples, respectively. Using the “SVA” software, we eliminated batch effects and normalized the data (17). Specifically, we re-analyzed the three datasets by combining the expression matrices from all three arrays before applying the normalization process. were then identified using the “Limma” package (18). The criteria for DEGs identification included an adjusted P-value < 0.05 and |logFC| > 1. Expression heat maps were generated using the “pheatmap” package.

The WGCNA method (19) is employed to explore the relationships between gene sets and clinical traits by constructing gene co-expression networks. Using the “WGCNA” R package, gene networks were constructed and modularized in several stages. WGCNA was applied to all genes that passed quality control (QC), not just DEGs. Initially, samples were clustered to identify significant outliers. Automated procedures were then used to construct the co-expression networks. Modules were identified through hierarchical clustering combined with dynamic tree cutting. To correlate modules with clinical traits, module membership (MM) and gene significance (GS) were calculated. Hub modules were defined by the highest Pearson correlation of MM and a P-value ≤ 0.05. High connectivity within modules and their clinical relevance were determined by MM values > 0.8 and GS values > 0.2, respectively. The network was constructed as an unsigned network to capture both positive and negative correlations between genes. The primary goal of this method is to identify key modules and genes that are strongly associated with clinical traits, which could serve as potential biomarkers for further investigation.

Functional enrichment analysis was performed to clarify the potential functions of prospective targets. Gene Set Enrichment Analysis (GSEA) (20) was performed using the immunologic signature gene set (C7 gene sets) from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/) using the “clusterProfiler” package. Gene ontology (GO) is commonly used to assign functions to genes, particularly in terms of biological pathways (BP), molecular functions (MF), and cellular components (CC). KEGG (21) enrichment analysis offers insights into gene functions and associated high-level genomic information. The R packages “GOplot” and “clusterProfiler” were used to explore the GO functions of potential mRNAs and improve KEGG pathways, leading to a better comprehension of the carcinogenic implications of target genes (22).

Feature genes, crucial for diagnosing RA, were isolated from previously identified genes. SVM-RFE, a machine learning technique, was employed to train a subset of features (genes) across different categories (such as patient vs. healthy groups), identifying the most predictive genes for RA diagnosis (23). The “glmnet” package in R was used to perform LASSO regression, which aims to select linear models by shrinking coefficients of less important variables to zero, thereby retaining the most relevant features for predictive modeling. LASSO classification, utilizing binomial distribution variables, was performed with a lambda value set at the 1-SE criterion, achieving optimal performance with only 10 cross-validation variables. RandomForest ranked the genes, with values above 0.25 considered significant (24). To further strengthen the analysis, advanced machine learning models, including CNN, were also applied. The CNN model was implemented using Python to capture complex patterns in gene expression data, offering a deep learning approach. The performance of CNN was compared to traditional machine learning models, assessing its ability to distinguish RA patients from healthy controls. The intersection of results from LASSO regression, SVM-RFE, RandomForest, and CNN analysis revealed the most significant feature genes in this study.

The “pROC” package was used to compute the area under the curve (AUC), generate Receiver Operating Characteristic (ROC) curves, and assess the diagnostic usefulness of signature genes in separating RA cases from normal samples using the GSE77298 dataset (25).

A nomogram model predicting RA recurrence was constructed using the “rms” package. In this model, “score” denotes the score for each variable, while “total score” represents the cumulative sum of all individual scores. Using calibration curves, the prediction accuracy of the model was evaluated. Clinical impact curves and decision curve analysis were used to assess the model’s clinical usefulness.

GeneMANIA (26) (http://www.genemania.org) enables the creation of Protein–Protein Interaction (PPI) networks, facilitating the prediction of gene functions and identification of genes with similar effects. Many bioinformatics techniques, including physical interactions, co-expression, co-localization, gene enrichment analysis, genetic relationships, and site predictions, are used by this network integration program. GeneMANIA was used in this research to examine the PPI networks including signature genes.

Single-sample gene set enrichment analysis (ssGSEA) (20) was used to quantify the infiltration scores of 16 immune types and the activities of 13 immune-related pathways. The infiltration score represents the relative abundance of each immune cell type or pathway in a given sample. The Spearman rank correlation coefficient was computed using the “corrplot” software to investigate the association between immunological state and certain genes.

Synovial tissue samples were collected from nine patients diagnosed with RA at Nanchong Central Hospital, all of whom met the clinical and imaging diagnostic criteria for RA. Additionally, nine young patients with no history of RA and undergoing surgical treatment for meniscal injuries, provided synovial tissue samples as trauma controls (TC). The inclusion criteria for the RA group required patients to meet both clinical and imaging diagnostic criteria for RA, while those with other autoimmune diseases or infections were excluded. For the TC group, inclusion was restricted to patients without a history of RA, and individuals with inflammatory joint diseases or autoimmune disorders were excluded. Sample IDs and demographic details for both groups are outlined in Supplementary Table 1.

To extract total RNA, synovial tissue was treated with Trizol reagent (Invitrogen, Carlsbad, CA, USA). After the RNA was separated, it was reverse transcribed into complementary DNA (cDNA) using gDNA Eraser and the PrimeScript RT Reagent Kit from Takara Bio, Inc. in Kyoto, Japan. SYBR Green qPCR Mix (Takara Bio, Inc., Japan) was used for qRT-PCR tests, which were performed using CFX96 quantitative PCR equipment (Bole, Inc., California, USA). Target gene expression levels were measured using the 2^-ΔΔCq method, normalized to GAPDH (27). Table 1 contains a collection of qRT-PCR primer sequences.

Proteins were extracted from synovial tissue using RIPA buffer (Beyotime, China) and quantified with BCA Protein Assay Kits (Beyotime). After that, the proteins were heated for 10 minutes and mixed in a 1:3 ratio with loading buffer to denature them. The proteins were separated using SDS-PAGE and then placed onto PVDF membranes (Millipore, USA). The PVDF membranes were subjected to an hour of blocking with 5% BSA at 4 °C. Subsequently, they were incubated with the primary antibodies against PCDH9 (Abcam, ab233710), GABARAPL1 (Abcam, ab109364), FKBP5 (Proteintech, 14155-1-AP), SLAMF8 (Novus Biologicals, AF1907), and GAPDH (Proteintech, 10494-1-AP) for an overnight period at 4°C. Membranes were incubated with the primary antibody for one hour at room temperature, and then treated with horseradish peroxidase-linked secondary antibody (Proteintech, 1:5000). An ECL kit (Biosharp, China) was used to observe the signal bands, and ImageJ software (Bethesda, MD, USA) was used to analyze the signal bands for gray values.

Synovial tissue samples were taken and preserved for 48 hours in 10% neutral buffered formalin. The paraffin-embedded fixed tissues were sectioned after being fixed. Hematoxylin and eosin were used to stain the sections after they had been deparaffinized in xylene and progressively dried in ethanol. Sections were counterstained with eosin after differentiation. After being dried in graded ethanol, the dyed slices were mounted using neutral resin. The specimens were then studied using an Olympus fluorescent microscope (Japan).

This investigation validated diagnostic biomarkers using immunofluorescence (IF) labeling. To retrieve antigens, deparaffinized synovial tissue slices were soaked in a 10 mM citrate buffer solution (Solarbio, China). The sections were incubated with the primary antibody solution overnight at 4°C after being blocked for an hour with 5% BSA (Solarbio). The slices were incubated with a fluorescence-labeled secondary antibody solution (Proteintech) for an extra hour the following day. After incubating the sections, they were cleaned three times for ten minutes each using 1× TBST solution. After the nuclei were counterstained with the blue fluorescent dye DAPI, the coverslip was sealed with an anti-fluorescence quencher. After that, the samples were examined and captured on camera using an Olympus fluorescent microscope (Japan).

To assess the clinical relevance of the identified biomarkers, we examined the associations between gene expression levels measured by quantitative PCR (qPCR) and serum inflammatory indices within the rheumatoid arthritis (RA) cohort. Synovial tissue–derived qPCR expression values for the candidate genes (e.g., GABARAPL1, FKBP5, PCDH9, and SLAMF8) were matched to the corresponding clinical data from the same patients. ESR (mm/h) and CRP (mg/L) values were collected from routine laboratory testing at the time of sample acquisition. Only RA patients with available ESR and CRP measurements were included in the correlation analysis; cases with missing values were excluded listwise.

Protein-drug interaction networks play a critical role in drug development (28). In this study, we used the Enrichr tool (29) to analyze transcriptomic signatures from the DSigDB database (30), identifying potential drug candidates targeting core genes. The ten most promising drugs were selected based on their P-values. These drugs were recommended for targeting hub genes due to their potential therapeutic benefits for RA.

Statistical analyses were performed using Prism 8.0 (GraphPad Software, San Diego, CA, USA). Normality was assessed using the Shapiro–Wilk test. For comparisons between two groups, an unpaired two-tailed Student’s t test was applied when data were normally distributed; otherwise, the non-parametric Mann–Whitney U test was used. For comparisons among multiple groups, one-way ANOVA with appropriate post hoc testing was performed for normally distributed data with homogeneity of variance; otherwise, the Kruskal–Wallis test followed by Dunn’s multiple-comparison test was applied. Data are presented as mean ± standard deviation (SD) for normally distributed variables and as median (interquartile range, IQR) for non-normally distributed variables. A two-sided p value < 0.05 was considered statistically significant. To ensure robustness and repeatability, all in vitro experiments were conducted in triplicate.

The research workflow is illustrated in Figure 1. Data standardization is shown in a box plot, where various colors represent distinct data sets, with rows corresponding to samples and columns indicating gene expression values (Figure 2A). Figure 2B demonstrates the PCA results from multiple data sets before addressing batch effects, using different colors to differentiate the data sets. Figure 2C shows the PCA results after batch effect removal, indicating that the integration of the three data sets allows for collective analysis. Applying the standards of |logFC| > 1 and an adjusted P-value < 0.05, 543 genes were found to be DEGs, of which 220 were down-regulated and 323 were up-regulated. Figure 2D shows the volcano plots of DEGs, and Figure 2E shows a heat map of the top 50 genes.

Datasets GSE12021, GSE55457, and GSE55235 were obtained from the GEO database, consisting of 29 normal and 35 RA samples. These samples were clustered using WGCNA to remove obvious outliers, with a defined threshold (Figure 3A). Figure 3B shows a soft-threshold power of 8, with an R² > 0.9 and strong average connectedness. A clustering height limit of 0.25 was set to merge strongly related modules, resulting in 21 modules for further analysis (Figure 3C). These modules were displayed on the clustering tree (Figure 3D). Module correlation analysis showed no significant associations among them (Figure 3E). The reliability of the module delineation was validated by transcription correlation analysis within the modules, which also revealed no significant connections (Figure 3F). To explore the relationship between modules and clinical features, correlations between module eigengene (ME) values and clinical features were analyzed. The pink module showed a positive correlation with normal samples (r = 0.83, P = 4e−17) and a negative correlation with RA samples (r = -0.83, P = 4e−17) (Figure 3G). This module, strongly associated with RA, was identified as clinically relevant, as shown in the MM vs. GS scatter plot (Figure 3H). Further analysis was then performed on all genes within the pink module.

A Venn diagram was used to overlap critical module genes with DEGs, resulting in the identification of 273 overlapping genes (Figure 4A). Functional analysis elucidated the biological functions of these DEGs within the modules. GSEA results indicated that these DEGs are associated with allograft rejection, autoimmune thyroid disease, pyruvate metabolism, and tyrosine metabolism (Figure 4B). KEGG analysis linked these DEGs to cytokine-cytokine receptor interactions, Th17 cell differentiation, chemokine signaling, and T cell receptor signaling pathways (Figure 4C). DO analysis demonstrated involvement of these DEGs in lymphoblastic leukemia, hepatitis, lung diseases, and chronic lymphocytic leukemia (Figure 4D). GO enrichment analysis revealed that the module DEGs enhance lymphocyte activation, mononuclear cell differentiation, leukocyte activation, T cell activation, and activities of chemokines, cytokines, and G protein-coupled receptors (Figure 4E).

To identify feature genes, we utilized key module genes from WGCNA analysis and 273 core genes derived from the intersection of DEGs as inputs for four advanced machine learning algorithms: LASSO regression, RandomForest, SVM-RFE, and CNN. These algorithms classified genes based on their ability to distinguish RA samples from healthy controls. The results of the SVM-RFE analysis are shown in Figures 5A, B, with the associated gene list provided in Supplementary Table 2. LASSO regression identified twelve genes based on statistically significant univariate factors, as shown in Figure 5C and Supplementary Table 3. RandomForest was used to assess the error rate, the number of classification trees, and the relevance of the 273 genes, incorporating feature selection. The results of this analysis are shown in Figures 5D, E with the relevant gene list available in Supplementary Table 4. The Performance Metrics of Machine Learning Models for RA Biomarker Prediction are included in Supplementary Table 5. To complement these methods, CNN was applied for further analysis. CNN captured complex gene expression patterns and improved classification accuracy. The CNN results, including the most significant feature genes, are presented in Figure 5F and detailed in Supplementary Table 6. Finally, a Venn diagram (Figure 5G) was used to identify eleven overlapping genes shared by all four methods, highlighting the most consistently identified feature genes.

To evaluate the diagnostic efficacy of EBF2, FKBP5, GABARAPL1, MREG, NEFM, PCDH9, RTN1, SFRP1, SLAMF8, SLC16A7, and SNX10, ROC analysis was conducted. The AUC values obtained were: EBF2 (0.795), FKBP5 (0.821), GABARAPL1 (0.955), MREG (0.848), NEFM (0.598), PCDH9 (0.857), RTN1 (0.688), SFRP1 (0.705), SLAMF8 (0.973), SLC16A7 (0.616), and SNX10 (0.768) (Figure 6). The five genes with the highest AUC values—FKBP5, GABARAPL1, MREG, PCDH9, and SLAMF8—were used to develop nomogram models for RA diagnostics using the “rms” package (Figure 7A). Their predictive accuracy was evaluated using calibration curves. Calibration curves showed minimal differences between actual and predicted RA risks, indicating high accuracy of the nomogram models (Figure 7B). The model’s accuracy was further verified through additional ROC curve analysis (Figure 7C). Validation with datasets GSE12021, GSE55235, and GSE55457 corroborated these findings (Figure 7D). These findings suggest that the five key genes are involved in the pathogenesis of RA.

Gene correlations were examined, revealing both positive and negative correlations among the expressions of FKBP5, GABARAPL1, MREG, PCDH9, and SLAMF8, as depicted in Figure 8A. This indicates significant functional relationships among these five genes. Subsequently, using the online platform GeneMANIA (http://genemania.org/), an intuitive network diagram was created to illustrate the interactions and relationships among these genes, highlighting their closely linked functions (Figure 8B).

Using ssGSEA, the relationship between immune infiltration in RA patients and healthy controls was investigated further. When statistically non-significant data were removed from the analysis, it was shown that although immune cell infiltration in mast cells was lower in RA patients than in controls, other immune-related pathway activities were higher in the overall RA group (Figure 9A). Subsequently, the “corrplot” package was utilized to analyze the correlation between signature genes and immune infiltration. FKBP5, GABARAPL1, and PCDH9 exhibited strong negative correlations with several immune functions, including CD8+ T cells, checkpoints, cytolytic activity, inflammation promotion, MHC class I, neutrophils, T-cell co-stimulation, Tfh, Th1/2 cells, TIL, Treg, and Type I IFN response. Conversely, MREG and SLAMF8 demonstrated strong positive correlations with functions such as Type I IFN response, Tfh, TIL, Th1/2 cells, T-cell co-inhibition and co-stimulation, checkpoints, and cytolytic activity (Figure 9B). These results imply that the hallmark genes may modify immunological mechanisms along the course of RA.

To identify potential biomarkers for RA, we selected five genes—FKBP5, GABARAPL1, MREG, PCDH9, and SLAMF8—based on their strong diagnostic potential and involvement in immune regulation, as indicated by previous bioinformatics analyses, suggesting their role in the disease’s development. Synovial tissues from nine patients diagnosed with RA and nine trauma control (TC) subjects were collected. Magnetic resonance imaging (MRI) analysis revealed synovial thickening in RA patients compared to the TC group (Figure 10A). H&E staining of RA synovial tissues indicated enhanced cellular density and heterogeneity in the intercellular matrix. Notably, significant increases in both the number and staining intensity of cell nuclei were observed, indicating elevated cell proliferation and substantial inflammatory cell infiltration (Figure 10B). qRT-PCR analyses demonstrated decreased expression of FKBH5, GABARAPL1, and PCDH9, and increased expression of SLAMF8 in the RA synovium; the rise in MREG expression was not statistically significant (Figure 10C). Western blot and IF analyses confirmed the downregulation of FKBH5, GABARAPL1, and PCDH9, along with the upregulation of SLAMF8, in the RA group compared to the TC group (Figures 10D, E). Spearman’s rank correlation analysis showed that FKBP5, GABARAPL1, and PCDH9 expression levels were negatively correlated with inflammatory activity, whereas SLAMF8 exhibited a positive correlation (Figure 11). All of these results show that the expression levels of FKBH5, GABARAPL1, PCDH9, and SLAMF8 in RA synovium vary significantly from those in TC synovium.

Understanding the structural factors influencing receptor sensitivity requires examining protein–drug interactions (31, 32). In this study, we used transcriptomic data from the DSigDB database and the Enrichr tool to identify 244 candidate pharmacological compounds associated with hub genes that may be relevant to RA. While these compounds have shown promise in in silico analyses as candidates for modulating hub gene activity, further validation and experimental studies are necessary to assess their potential role in RA treatment. The top ten compounds identified from the DSigDB database, ranked by P-values, are listed in Supplementary Table 7. It is important to note that many of these compounds, such as daunorubicin, cupric oxide, and methaneseleninic acid, are cytotoxic or have been studied in contexts other than RA. Therefore, the current evidence supporting their relevance to RA is computational and should be considered hypothesis-generating.