Abstract
Objective:
Rheumatoid arthritis (RA) is a chronic autoimmune disorder that significantly impacts quality of life. Despite extensive research, its pathogenesis remains unclear. This study aims to identify potential diagnostic biomarkers and therapeutic targets for RA.
Methods:
This study integrated patient data from three Gene Expression Omnibus (GEO) databases to analyze gene expression in RA. Using Weighted Gene Correlation Network Analysis (WGCNA), we identified key genes, which were then compared with differentially expressed genes (DEGs) to uncover RA-related genes. Functional enrichment analysis provided insights into the biological roles of these genes. To refine our findings, we applied three algorithms—RandomForest, SVM-REF, LASSO, and Convolutional Neural Networks (CNN)—to pinpoint a subset of core genes. We evaluated their diagnostic potential through receiver operating characteristic (ROC) curves and selected the top five genes with the highest area under the curve (AUC) values for constructing a predictive nomogram model. An interaction analysis was performed to investigate the relationship between these genes and immune cell infiltration. Finally, the expression of these core genes was validated in the synovial tissues of RA patients. Drug-protein interaction relationships were predicted using the DSigDB database.
Results:
Differential expression analysis identified 543 DEGs. We subsequently applied WGCNA to compare these DEGs with significant module genes, resulting in the identification of 273 key genes. Functional enrichment analysis indicated that these genes were primarily involved in inflammatory response pathways. Further analysis using four machine learning algorithms identified 11 core genes. Of these, the five genes with the highest AUC values were selected to construct a robust nomogram model. Immune infiltration analysis revealed significant differences in immune cell levels and pathways between RA patients and healthy controls, which were correlated with the expression of these five genes. Validation through quantitative real-time PCR (qRT-PCR), Western blot, and immunofluorescence (IF) confirmed that GABARAPL1, FKBP5, and PCDH9 expression was lower in RA synovial tissues compared to healthy controls, while SLAMF8 expression was elevated. Additionally, potential therapeutic drugs targeting these key genes, including (+)-chelidonine, daunorubicin, and bisacodyl, were predicted.
Conclusion:
GABARAPL1, FKBP5, PCDH9, and SLAMF8 are identified as potential biomarkers for RA, offering insights into future therapeutic strategies.
Introduction
Rheumatoid arthritis (RA) is a chronic autoimmune disease marked by systemic inflammation and persistent synovitis (1–3). The pathology of RA involves immune cell infiltration, hyperplasia of the synovial lining, pannus formation, and destruction of articular cartilage and bone (4). Currently, RA affects approximately 0.5–1% of adults worldwide, imposing significant economic burdens on both society and patients (5). Females are particularly more susceptible to RA, with a risk two to three times higher than that of males (6). Symptoms of RA include morning stiffness, which, if untreated, can progress to focal necrosis, granulation adhesion, and fibrous tissue formation on joint surfaces. This progression may lead to severe joint ankylosis, destruction, deformities, and disability (7). Researches show genetic, environmental, and immune factors are recognized to play significant roles (1, 8), despite extensive research efforts, the complete pathogenesis of RA remains elusive. Furthermore, current diagnostic approaches for RA are inadequate, underscoring the critical need for the development of reliable and efficient biomarkers to improve patient outcomes (9).
The etiology and pathogenesis of RA involve a multifaceted interplay of infectious agents, genetic susceptibilities, and environmental factors that trigger immune responses, resulting in synovitis. Studies indicate that the progression of RA is greatly affected by the atypical morphology and gene expression observed in RA macrophages (RASM) and synovial fibroblasts (RASF) (10). B cells exacerbate RA by producing proteins such as rheumatoid factor (RF), anti-citrullinated protein antibodies (ACPA), and pro-inflammatory cytokines. These proteins form immune complexes with self-antigens, thereby sustaining the disease process (11). Upon stimulation by antigen-presenting cells, CD4+ T cells evolve into five distinct lineages of T helper (Th) cells—Th17, Th2, Th1, T follicular helper (Tfh), and T regulatory (Treg) cells (12). In RA, these T cells predominantly activate macrophages and fibroblasts, promoting their transformation into cells that damage tissues (13). Despite this understanding, the precise mechanisms underlying gene and protein expression in RA-associated synovium remain poorly understood.
The primary therapeutic options for RA include anti-inflammatory drugs, analgesics, and disease-modifying antirheumatic drugs (DMARDs). Anti-inflammatory drugs and analgesics only relieve RA symptoms without stopping disease progression. Because immune responses are crucial in the pathogenesis of RA, DMARDs are the preferred treatment method. Although DMARDs effectively reduce RA symptoms, many patients still experience treatment failures, including nonresponse or limited effectiveness (14). Over the last decade, biological agents have increasingly been tested in clinical trials. These agents target immune cells for immunomodulation and are often used alongside DMARDs to manage RA (15). Whether used alone or combined with new biologics, the optimal treatment regimen for RA is still being researched. Thus, exploring genes linked to the onset and progression of RA is essential for enhancing early diagnosis and developing better therapeutic strategies.
Substantial attention is currently focused on identifying disease-specific biomarkers using bioinformatics and genome sequencing technologies (16). Our research aimed to pinpoint key biomarkers related to DEGs and immune infiltration in RA, with the ultimate goal of developing novel diagnostic and therapeutic strategies for this condition.
Materials and methods
Data processing
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.
Construction of WGCNA
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 analysis
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).
Identification of hub genes
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.
Curve analysis of ROC
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).
Development and validation of diagnostic nomogram
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.
Construction of PPI network
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.
Immune infiltration analysis using ssGSEA
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.
Clinical specimens
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.
RNA extraction and qRT-PCR
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.
Table 1
| Gene | Forward and reverse primer |
|---|---|
| FKBP5 | F:5’ GCGAAGGAGAAGACCACGACAT 3’ |
| R:5’ TAGGCTTCCCTGCCTCTCCAAA 3’ | |
| GABARAPL1 | F:5’ CCCTCCCTTGTTATCATCCA 3’ |
| R:5’ ACTCCCACCCCACAAAATCC 3’ | |
| MREG | F:5’ CCCTTGGCATTTTATCTGGA 3’ |
| R:5’ AAGCTGCATTCACAGCATTG 3’ | |
| PCDH9 | F:5’ TCCCAACTCTGATGGGCCTTTGGG 3’ |
| R:5’ GGCTCTGGTCAGGGTGTGCC 3’ | |
| SLAMF8 | F:5’ CTGATGGTGGATACAAGGG 3’ |
| R:5’ GGAAATGGACGTAACGGA 3’ | |
| GAPDH | F:5’ TGTGTCCGTCGTGGATCTGA 3’ R:5’ TTGCTGTTGAAGTCGCAGGAG 3’ |
The sequence of primers.
Protein extraction and Western blot
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.
Hematoxylin-eosin staining
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).
Immunofluorescence staining
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).
Clinical correlation analysis
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.
Evaluation of applicant drugs
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 analysis
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.
Results
DEGs screening and data preprocessing
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.
Figure 1
Figure 2
Data collection and module analysis
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 R2 > 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.
Figure 3
Functional analysis of critical module genes and DEGs
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).
Figure 4
Selection of feature genes
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.
Figure 5
Modeling and testing of an RA diagnostic nomogram
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.
Figure 6
Figure 7
Trait gene interaction analysis
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).
Figure 8
ssGSEA analysis and correlation analysis
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.
Figure 9
Validation of the hub genes expression
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.
Figure 10
Figure 11
Identification of candidate drugs
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.
Discussion
Rheumatoid arthritis (RA) severely impairs the quality of life for sufferers and causes irreparable damage to the synovial lining of joints (33). The incidence of RA has grown during the last several decades despite advancements in therapy techniques (34). Early diagnosis and treatment of RA are crucial for preventing joint damage and improving quality of life. However, the lack of effective biomarkers makes diagnosing early-stage RA challenging. Identifying novel, reliable biomarkers is essential for RA therapy and early diagnosis. The purpose of this research was to find new biomarkers for RA diagnosis and treatment and to look at how they relate to immune infiltration. Several bioinformatics techniques and in vivo tests were used to identify and verify four diagnostic potential diagnostic biomarkers for RA. These four diagnostic potential diagnostic biomarkers showed notable variations in expression within the immune microenvironment.
In this study, we identified 543 DEGs by comparing gene expression in samples from the healthy controls and RA groups. We then performed WGCNA analysis, which revealed 21 co-expression modules. Genes from modules highly correlated with RA phenotypes were selected and intersected with the DEGs, resulting in 273 key genes. A subsequent GO enrichment analysis indicated that these key genes were predominantly associated with cytokine receptor function and the activation of inflammatory cells. Additionally, KEGG enrichment analysis demonstrated their involvement in T-cell receptor signaling, chemokine signaling, and cytokine-cytokine receptor interactions. The enrichment analysis further revealed that RA synovium exhibits strong immune activation and signal transduction effects, which play major roles in RA-related synovial inflammation, arthritis, and joint pain. Arthritis and joint pain are well-established as the primary clinical manifestations of RA (35).
Machine learning provides an objective approach for predicting patient status, simultaneously revealing previously unknown interactions and identifying novel biological signatures (36, 37). Using four machine learning algorithms and ROC curves, we identified GABARAPL1, FKBP5, MREG, PCDH9, and SLAMF8 as diagnostic signature biomarkers. As demonstrated by the diagnostic nomogram, these five RA biomarkers possess significant diagnostic value. Furthermore, clinical tissue validation data showed notable differences in the expression of four diagnostic biomarkers (GABARAPL1, FKBP5, PCDH9, and SLAMF8) between TC and RA tissues. Immune infiltration analysis also revealed strong correlations between these biomarkers and various immune cell functions, as well as immune-related pathways. Finally, we identified potential drugs targeting these four DEGs, including (+)-chelidonine, daunorubicin, and bisacodyl. Investigating the interactions between these key genes and their respective drugs could significantly contribute to advancing drug development for RA. Therefore, based on our findings, we propose that GABARAPL1, FKBP5, PCDH9, and SLAMF8 may serve as biomarkers for the early diagnosis of RA.
GABARAPL1 is a member of the GABARAP family and is also known as autophagy-related 8 (ATG8) or glandular epithelial cell protein 1 (GEC1) (38). In addition to being engaged in the movement of proteins and vesicles, the GABARAP group members (GABARAP, GABARAPL1, and GABARAPL2) and their homologs (LC3 and Atg8) are also involved in many other processes, including autophagy, apoptosis, cell development, and cancer (38). GEC1 regulates mitochondrial homeostasis, which is necessary for autophagosome maturation, as An et al. (39) showed. Further evidence that GEC1 maintains appropriate autophagic flow and is essential for controlling cell bioenergetics and metabolism was provided by Boyer-Guittaut et al. (40). GABARAPL1 may affect cellular energy and cytoskeletal reorganizations required for immune cell migration under its function in autophagy (41). Research suggests that GABARAPL1 may indirectly influence the threshold for immune cell death by regulating the autophagic response. For instance, increased autophagic activity mediated by GABARAPL1 might lessen the tendency for cell apoptosis and increase cell survival under stress (42, 43). According to our research, the RA synovium has downregulated GABARAPL1 expression. Furthermore, GABARAPL1’s ROC curve demonstrated that it has a strong diagnostic value for RA (AUC = 0.955). We believe that GABARAPL1 is a useful biomarker for the diagnosis of RA.
FKBP5, also known as FKBP51, is a member of the FK506 binding protein (FKBP) family and is involved in many biological functions, including transport, immunological control, and protein folding (44, 45). FKBP51 functions as a chaperone protein, affecting the production of proinflammatory cytokines by modulating the transcription factor nuclear factor kappa B (NF-κB) (46). Fascinatingly, Nakamura et al. (47) recently reported that RA patients’ bone marrow (BM) mononuclear cells had elevated expression levels of many genes, including FKBP5. Nonetheless, we discovered that the RA synovium had downregulated FKBP5 expression. Consequently, FKBP5 could be important for the development of RA illness.
PCDH9 belongs to the biggest subfamily in the cadherin group, the protocadherin family (48). A member of the δ1-subfamily that also contains PCDH1, PCDH7, PCDH9, and PCDH11, PCDH9 is involved in the formation and breakdown of cell adhesions (49). Protocadherins influence cell migration, survival, and growth inhibition, all of which are important aspects of carcinogenesis (50). Due to their excessive activation and tumor-like features, FLSs play a critical role in the onset and evolution of RA by producing matrix metalloproteinase (MMP), pro-inflammatory cytokine release, invasion, and aberrant proliferation (51–54). Thus, we speculate that PCDH9 plays a role in developing features resembling tumors in synoviocytes. Furthermore, we discovered that PCDH9 had a good diagnostic value for RA (AUC = 0.857) and was only weakly expressed in the synovium of RA patients. In our opinion, PCDH9 is a new and useful biomarker for diagnosing RA.
The signaling lymphocytic activation molecule family (SLAMF) is a family of surface receptors that regulate the activity of immune cells and is extensively expressed in hematopoietic cells, including T cells, NK cells, and macrophages (55). It has been shown that SLAMF8, a member of the SLAMF family, lowers ROS generation and inhibits macrophage activity (56). Additionally, it has been shown to be a noteworthy indicator of T-cell infiltration into malignancies (57). The TLR4/NF-κB signaling pathway was blocked by targeted suppression of SLAMF8, which decreased the severity of RA (58). Furthermore, we discovered that SLAMF8 had a strong diagnostic value for RA (AUC = 0.973) and was substantially expressed in the synovium of RA patients. Therefore, we think SLAMF8 is a useful biomarker for RA diagnosis.
Limitations
This research has certain limitations, even though our nomogram model, which is based on five distinctive genes, shows strong diagnostic predictive capacity for RA patients. First of all, certain demographic and clinical characteristics (e.g., sex, age, comorbidities) were not systematically evaluated in this database-derived study, highlighting the need for further clinical investigations to validate these findings. Secondly, the biological functions of the genes that have been discovered and their associations with RA remain incompletely understood. Lastly, biased validation data arise from the causes of clinical sample limits. To further confirm the practical applicability value of our nomogram model in identifying RA patients, we want to increase clinical patient recruitment.
Conclusion
In conclusion, we found that FKBP5, GABARAPL1, PCDH9, and SLAMF8 are diagnostic biomarkers for RA linked to immune infiltration via our investigation using a machine learning method. These results provide a fresh viewpoint on the pathophysiology of RA.
Statements
Data availability statement
The NCBI GEO repository has the datasets that this study is based on: GSE12021, GSE55235, GSE55457, GSE77298, and https://www.ncbi.nlm.nih.gov/geo/. The corresponding author may provide the experimental data utilized in this work upon reasonable request.
Ethics statement
The studies involving humans were approved by the Ethics Committee of Nanchong Central Hospital. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.
Author contributions
Z-WF: Writing – review & editing, Methodology, Investigation, Visualization, Software, Project administration, Funding acquisition, Validation, Resources, Writing – original draft, Data curation, Formal Analysis, Supervision, Conceptualization. M-KY: Writing – review & editing, Data curation, Methodology, Formal analysis, Project administration. X-DJ: Writing – review & editing, Data curation, Methodology, Formal analysis, Project administration. FY: Writing – review & editing, Writing – original draft. M-GG: Writing – review & editing, Writing – original draft. FC: Writing – original draft, Writing – review & editing. WL: Data curation, Software, Conceptualization, Visualization, Methodology, Investigation, Validation, Funding acquisition, Writing – review & editing, Supervision, Project administration, Resources, Formal Analysis, Writing – original draft. C-FY: Investigation, Project administration, Resources, Conceptualization, Visualization, Validation, Funding acquisition, Methodology, Supervision, Formal Analysis, Software, Writing – review & editing, Writing – original draft, Data curation.
Funding
The author(s) declared that financial support was received for this work and/or its publication. Northern Sichuan Health Humanities Research Special Project (NC25CB64); Special Scientific Research Project of Orthopedics (Shangantong) of Sichuan Medical Association (2023SAT30); Northern Sichuan Medical College-level Scientific Research Development Plan Project (CBY23-QNB04); Science and Technology Project of Sichuan Provincial Health Commission (24WSXT088); Project of Sichuan Provincial Medical Science and Technology Innovation Research Association (YCH-KY-YCZD2024-078); Real-world research project on the clinical diagnosis and treatment of spinal degenerative diseases (WKZX2024JZ0218).
Acknowledgments
This work is equally credited to Z-WF and M-KY as principal writers.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2026.1645257/full#supplementary-material
References
1
FiresteinGSMcInnesIB. Immunopathogenesis of rheumatoid arthritis. Immunity. (2017) 46:183–96. doi: 10.1016/j.immuni.2017.02.006
2
FengZWYangCFXiaoHFYuanFChenFZhangBet al. YTHDC1 regulates the migration, invasion, proliferation, and apoptosis of rheumatoid fibroblast-like synoviocytes. Front Immunol. (2024) 15:1440398. doi: 10.3389/fimmu.2024.1440398
3
FengZWTangYCShengXYWangSHWangYBLiuZCet al. Screening and identification of potential hub genes and immune cell infiltration in the synovial tissue of rheumatoid arthritis by bioinformatic approach. Heliyon. (2023) 9:e12799. doi: 10.1016/j.heliyon.2023.e12799
4
SmolenJSAletahaDBartonABurmesterGREmeryPFiresteinGSet al. Strand V et al: Rheumatoid arthritis. Nat Rev Dis Primers. (2018) 4:18001. doi: 10.1038/nrdp.2018.1
5
ScottDLWolfeFHuizingaTW. Rheumatoid arthritis. Lancet. (2010) 376:1094–108. doi: 10.1016/S0140-6736(10)60826-4
6
NgoSTSteynFJMcCombePA. Gender differences in autoimmune disease. Front Neuroendocrinol. (2014) 35:347–69. doi: 10.1016/j.yfrne.2014.04.004
7
TuressonCO’FallonWMCrowsonCSGabrielSEMattesonEL. Extra-articular disease manifestations in rheumatoid arthritis: incidence trends and risk factors over 46 years. Ann Rheum Dis. (2003) 62:722–7. doi: 10.1136/ard.62.8.722
8
JangSKwonEJLeeJJ. Rheumatoid arthritis: pathogenic roles of diverse immune cells. Int J Mol Sci. (2022) 23. doi: 10.3390/ijms23020905
9
AbdelhafizDBakerTGlascowDAAbdelhafizA. Biomarkers for the diagnosis and treatment of rheumatoid arthritis - a systematic review. Postgrad Med. (2023) 135:214–23. doi: 10.1080/00325481.2022.2052626
10
ElshabrawyHAChenZVolinMVRavellaSVirupannavarSShahraraS. The pathogenic role of angiogenesis in rheumatoid arthritis. Angiogenesis. (2015) 18:433–48. doi: 10.1007/s10456-015-9477-2
11
YapHYTeeSZWongMMChowSKPehSCTeowSY. Pathogenic role of immune cells in rheumatoid arthritis: implications in clinical treatment and biomarker development. Cells. (2018) 7. doi: 10.3390/cells7100161
12
Suárez-FueyoABradleySJTsokosGC. T cells in systemic lupus erythematosus. Curr Opin Immunol. (2016) 43:32–8. doi: 10.1016/j.coi.2016.09.001
13
ChoyE. Understanding the dynamics: pathways involved in the pathogenesis of rheumatoid arthritis. Rheumatol (Oxford). (2012) 51 Suppl 5:v3–11. doi: 10.1093/rheumatology/kes113
14
MittalNMittalRSharmaAJoseVWanchuASinghS. Treatment failure with disease-modifying antirheumatic drugs in rheumatoid arthritis patients. Singapore Med J. (2012) 53:532–6.
15
AbbasiMMousaviMJJamalzehiSAlimohammadiRBezvanMHMohammadiHet al. Strategies toward rheumatoid arthritis therapy; the old and the new. J Cell Physiol. (2019) 234:10018–31. doi: 10.1002/jcp.27860
16
LombeCPMeyerMPretoriusA. Bioinformatics prediction and analysis of microRNAs and their targets as biomarkers for prostate cancer: A preliminary study. Mol Biotechnol. (2022) 64:401–12. doi: 10.1007/s12033-021-00414-8
17
LeekJTJohnsonWEParkerHSJaffeAEStoreyJD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. (2012) 28:882–3. doi: 10.1093/bioinformatics/bts034
18
RitchieMEPhipsonBWuDHuYLawCWShiWet al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47–7. doi: 10.1093/nar/gkv007
19
LangfelderPHorvathS. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. (2008) 9:559. doi: 10.1186/1471-2105-9-559
20
SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMAet al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102
21
KanehisaMFurumichiMSatoYMatsuuraYIshiguro-WatanabeM. KEGG: biological systems database as a model of the real world. Nucleic Acids Res. (2025) 53:D672–d677. doi: 10.1093/nar/gkae909
22
WalterWSánchez-CaboFRicoteM. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics. (2015) 31:2912–4. doi: 10.1093/bioinformatics/btv300
23
SanzHValimCVegasEOllerJMReverterF. SVM-RFE: selection and visualization of the most relevant features through non-linear kernels. BMC Bioinf. (2018) 19:432. doi: 10.1186/s12859-018-2451-4
24
IshwaranHKogalurUB. Consistency of random survival forests. Stat Probab Lett. (2010) 80:1056–64. doi: 10.1016/j.spl.2010.02.020
25
RobinXTurckNHainardATibertiNLisacekFSanchezJCet al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinf. (2011) 12:77. doi: 10.1186/1471-2105-12-77
26
FranzMRodriguezHLopesCZuberiKMontojoJBaderGDet al. GeneMANIA update 2018. Nucleic Acids Res. (2018) 46:W60–w64. doi: 10.1093/nar/gky311
27
LivakKJSchmittgenTD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262
28
YanWZhangDShenCLiangZHuG. Recent advances on the network models in target-based drug discovery. Curr Top Med Chem. (2018) 18:1031–43. doi: 10.2174/1568026618666180719152258
29
KuleshovMVJonesMRRouillardADFernandezNFDuanQWangZet al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. (2016) 44:W90–7. doi: 10.1093/nar/gkw377
30
YooMShinJKimJRyallKALeeKLeeSet al. DSigDB: drug signatures database for gene set analysis. Bioinformatics. (2015) 31:3069–71. doi: 10.1093/bioinformatics/btv313
31
Al-MustanjidMMahmudSHRoyelMRIRahmanMHIslamTRahmanMRet al. Detection of molecular signatures and pathways shared in inflammatory bowel disease and colorectal cancer: a bioinformatics and systems biology approach. Genomics. (2020) 112:3416–26. doi: 10.1016/j.ygeno.2020.06.001
32
MahmudSHChenWMengHJahanHLiuYHasanSM. Prediction of drug-target interaction based on protein features using undersampling and feature selection techniques with boosting. Analytical Biochem. (2020) 589:113507. doi: 10.1016/j.ab.2019.113507
33
ChowdhuryTDuttaJNoelPIslamRGonzalez-PeltierGAzadSet al. An overview on causes of nonadherence in the treatment of rheumatoid arthritis: its effect on mortality and ways to improve adherence. Cureus. (2022) 14:e24520. doi: 10.7759/cureus.24520
34
DuongSQCrowsonCSAthreyaAAtkinsonEJDavisJM3rdWarringtonKJet al. Clinical predictors of response to methotrexate in patients with rheumatoid arthritis: a machine learning approach using clinical trial data. Arthritis Res Ther. (2022) 24:162. doi: 10.1186/s13075-022-02851-5
35
SmolenJSAletahaDMcInnesIB. Rheumatoid arthritis. Lancet. (2016) 388:2023–38. doi: 10.1016/S0140-6736(16)30173-8
36
ErringtonNIremongerJPickworthJAKariotisSRhodesCJRothmanAMet al. A diagnostic miRNA signature for pulmonary arterial hypertension using a consensus machine learning approach. EBioMedicine. (2021) 69:103444. doi: 10.1016/j.ebiom.2021.103444
37
NeumannURiemenschneiderMSowaJPBaarsTKälschJCanbayAet al. Compensation of feature selection biases accompanied with improved predictive performance for binary classification by using a novel ensemble feature selection approach. BioData Min. (2016) 9:36. doi: 10.1186/s13040-016-0114-4
38
Le GrandJNChakramaFZSeguin-PySFraichardADelage-MourrouxRJouvenotMet al. GABARAPL1 (GEC1): original or copycat? Autophagy. (2011) 7:1098–107. doi: 10.4161/auto.7.10.15904
39
AnHKChungKMParkHHongJGimJEChoiHet al. CASP9 (caspase 9) is essential for autophagosome maturation through regulation of mitochondrial homeostasis. Autophagy. (2020) 16:1598–617. doi: 10.1080/15548627.2019.1695398
40
Boyer-GuittautMPoilletLLiangQBôle-RichardEOuyangXBenavidesGAet al. The role of GABARAPL1/GEC1 in autophagic flux and mitochondrial quality control in MDA-MB-436 breast cancer cells. Autophagy. (2014) 10:986–1003. doi: 10.4161/auto.28390
41
ChoiYBowmanJWJungJU. Autophagy during viral infection - a double-edged sword. Nat Rev Microbiol. (2018) 16:341–54. doi: 10.1038/s41579-018-0003-6
42
SahaSPanigrahiDPPatilSBhutiaSK. Autophagy in health and disease: A comprehensive review. BioMed Pharmacother. (2018) 104:485–95. doi: 10.1016/j.biopha.2018.05.007
43
DereticVSaitohTAkiraS. Autophagy in infection, inflammation and immunity. Nat Rev Immunol. (2013) 13:722–37. doi: 10.1038/nri3532
44
BaughmanGWiederrechtGJCampbellNFMartinMMBourgeoisS. FKBP51, a novel T-cell-specific immunophilin capable of calcineurin inhibition. Mol Cell Biol. (1995) 15:4395–402. doi: 10.1128/MCB.15.8.4395
45
MendonçaMSMangiavacchiPMRios ÁFL. Regulatory functions of FKBP5 intronic regions associated with psychiatric disorders. J Psychiatr Res. (2021) 143:1–8. doi: 10.1016/j.jpsychires.2021.08.014
46
ZannasASJiaMHafnerKBaumertJWiechmannTPapeJCet al. Epigenetic upregulation of FKBP5 by aging and stress contributes to NF-κB-driven inflammation and cardiovascular risk. Proc Natl Acad Sci U.S.A. (2019) 116:11370–9. doi: 10.1073/pnas.1816847116
47
NakamuraNShimaokaYTouganTOndaHOkuzakiDZhaoHet al. Isolation and expression profiling of genes upregulated in bone marrow-derived mononuclear cells of rheumatoid arthritis patients. DNA Res. (2006) 13:169–83. doi: 10.1093/dnares/dsl006
48
HodisEWatsonIRKryukovGVAroldSTImielinskiMTheurillatJPet al. : A landscape of driver mutations in melanoma. Cell. (2012) 150:251–63. doi: 10.1016/j.cell.2012.06.024
49
PeekSLMahKMWeinerJA. Regulation of neural circuit formation by protocadherins. Cell Mol Life Sci. (2017) 74:4133–57. doi: 10.1007/s00018-017-2572-3
50
ChenYXiangHZhangYWangJYuG. Loss of PCDH9 is associated with the differentiation of tumor cells and metastasis and predicts poor survival in gastric cancer. Clin Exp Metastasis. (2015) 32:417–28. doi: 10.1007/s10585-015-9712-7
51
NygaardGFiresteinGS. Restoring synovial homeostasis in rheumatoid arthritis by targeting fibroblast-like synoviocytes. Nat Rev Rheumatol. (2020) 16:316–33. doi: 10.1038/s41584-020-0413-5
52
KomatsuNTakayanagiH. Mechanisms of joint destruction in rheumatoid arthritis - immune cell-fibroblast-bone interactions. Nat Rev Rheumatol. (2022) 18:415–29. doi: 10.1038/s41584-022-00793-5
53
BartokBFiresteinGS. Fibroblast-like synoviocytes: key effector cells in rheumatoid arthritis. Immunol Rev. (2010) 233:233–55. doi: 10.1111/j.0105-2896.2009.00859.x
54
WuZMaDYangHGaoJZhangGXuKet al. Fibroblast-like synoviocytes in rheumatoid arthritis: Surface markers and phenotypes. Int Immunopharmacol. (2021) 93:107392. doi: 10.1016/j.intimp.2021.107392
55
ChenSCaiCLiZLiuGWangYBlonskaMet al. Dissection of SAP-dependent and SAP-independent SLAM family signaling in NKT cell development and humoral immunity. J Exp Med. (2017) 214:475–89. doi: 10.1084/jem.20161312
56
WangGAbadía-MolinaACBergerSBRomeroXO’KeeffeMSRojas-BarrosDIet al. Cutting edge: Slamf8 is a negative regulator of Nox2 activity in macrophages. J Immunol. (2012) 188:5829–32. doi: 10.4049/jimmunol.1102620
57
SugimotoAKataokaTRItoHKitamuraKSaitoNHirataMet al. SLAM family member 8 is expressed in and enhances the growth of anaplastic large cell lymphoma. Sci Rep. (2020) 10:2505. doi: 10.1038/s41598-020-59530-1
58
QinWRongXYuCJiaPYangJZhouG. Knockout of SLAMF8 attenuates collagen-induced rheumatoid arthritis in mice through inhibiting TLR4/NF-κB signaling pathway. Int Immunopharmacol. (2022) 107:108644. doi: 10.1016/j.intimp.2022.108644
Summary
Keywords
diagnostic biomarker, immune infiltration, machine learning, rheumatoid arthritis, WGCNA
Citation
Feng Z, Yang M, Jia X, Yuan F, Guo M, Chen F, Li W and Yang C (2026) Unveiling potential diagnostic biomarkers for rheumatoid arthritis through integrated gene expression analysis. Front. Immunol. 17:1645257. doi: 10.3389/fimmu.2026.1645257
Received
11 June 2025
Revised
05 January 2026
Accepted
09 January 2026
Published
24 February 2026
Volume
17 - 2026
Edited by
Pier Paolo Sainaghi, University of Eastern Piedmont, Italy
Reviewed by
Mu Kai-lang, Guizhou University of Traditional Chinese Medicine, China
Haseeb Nisar, King Fahd University of Petroleum and Minerals, Saudi Arabia
Updates
Copyright
© 2026 Feng, Yang, Jia, Yuan, Guo, Chen, Li and Yang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chen-fei Yang, 389969499@qq.com
†These authors have contributed equally to this work
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.