Machine learning-based identification of CXCL13 and CXCL10 as key biomarkers in rheumatoid arthritis

Rheumatoid arthritis (RA) is one of the autoimmune diseases that involve joint damage and synovial membrane inflammation. There have been many attempts to examine the utility of biomarkers in the diagnosis of RA, but there is no clinically accepted biomarker for the diagnosis of RA due to its heterogeneity. Thus, in this study, we tried to find RA biomarkers using machine learning and bioinformatics approaches. Three corresponding gene expression datasets (GSE55457, GSE12021, and GSE77298) were downloaded from the Gene Expression Omnibus database. Batch effect correction was carried out using the ComBat function of the sva package in R. Differentially expressed genes (DEGs) compared to controls were identified using limma. To identify the biological functions and pathways of DEGs, Gene Ontology, KEGG, and gene set enrichment analysis (GSEA) were conducted. The protein–protein interaction (PPI) network was constructed in Cytoscape using the STRING database. Genes identified through PPI analysis and further filtered by LASSO and SVM-RFE regressions were considered as diagnostic biomarkers. Additional analysis of hub genes was conducted using the GSE89408 dataset. CIBERSORTx was used to evaluate the immune cell infiltration in RA. A total of 440 DEGs were selected, including 197 downregulated and 243 upregulated genes. The DEGs were involved in biological processes and pathways associated with the immune system and were suggested by functional enrichment analysis. GSEA enriched in cytokine–cytokine receptor interaction, hematopoietic cell lineage, primary immunodeficiency, and RA pathways. The hub genes CXCL13 and CXCL10 were chosen with the help of a machine learning algorithm. An external independent validation of the hub genes was conducted. In addition, ROC analysis showed a good diagnostic value of the two genes. The analysis of immune cell infiltration showed CXCL13 was associated with multiple immune-cell populations, including plasma cells, T-cell subsets, and M1 macrophages, whereas CXCL10 showed significant associations with M1 macrophages, plasma cells, and selected T-cell subsets. Correlation analysis showed a high correlation between chemokines CXCL13 and CXCL10 and many infiltrated immune cells. The present study indicates the potential usefulness of CXCL13 and CXCL10 genes as diagnostic biomarkers in RA. Synovial inflammation in RA is associated with the coordinated action of many genes involved in the recruitment of various immune cells to the inflammatory site. Although the term chemokines were used to designate them, these genes are the genes encoding molecules associated with the migration of particular immune cell types, and CXCR5 and CXCR3 may be considered the potential targets for future research.

Rheumatoid arthritis is a chronic systemic autoimmune disease characterized by persistent synovial inflammation leading to both local joint damage and systemic manifestations (Firestein, 2003). Synovial hyperplasia in RA is associated with pannus formation, which invades cartilage and bone, ultimately leading to progressive and irreversible joint destruction (Paleolog, 2002). In addition, RA is associated with extra-articular complications affecting organs such as the lungs, heart, and vasculature (Conforti et al., 2021). RA represents one of the leading autoimmune disorders on a global scale, with population-based studies estimating that approximately 0.4%–1% of the global population is affected (Deshmukh, 2023).

The identification of early disease activity and timely intervention is needed to prevent irreversible structural damage in joints and disability. However, currently used biomarkers, including rheumatoid factor and anti-cyclic citrullinated peptide antibodies, demonstrate limited sensitivity and fail to reliably reflect fluctuations in disease activity. These limitations highlight the need for more reliable and robust molecular biomarkers for RA diagnosis and progression monitoring (Schellek et al., 2000; Combe, 2007).

Rapid development in genomic platforms and technologies, such as microarrays, has enabled comprehensive genome-wide expression analysis (Satam et al., 2023). The integration of publicly available repositories such as the Gene Expression Omnibus (GEO) enhances statistical power and reduces batch effects, facilitating robust identification of differentially expressed genes (DEGs) and disease-associated molecular signatures in RA (Edgar et al., 2002).

RA pathogenesis is strongly associated with immune dysregulation within the synovial microenvironment. It is often characterized by T lymphocytes, B lymphocytes, macrophages, and dendritic cells. These resident immune cells secrete a range of pro-inflammatory cytokines, including TNF-α, IL-1β, GM-CSF, and interferon-γ, along with chemokines such as CXCL9, CXCL10, and CXCL13, which are involved in immune cell recruitment and retention in inflamed synovial tissue (Firestein and McInnes, 2017). Here, we used computational approaches such as CIBERSORTx, which have been widely used to characterize immune subset composition from data patterns in RA tissues (Zhou et al., 2021).

In this article, we integrated three gene expression datasets and applied bioinformatics workflows and machine learning classification approaches to identify potential diagnostic biomarkers for RA. Furthermore, we aimed to elucidate the underlying disease mechanisms with a particular focus on immune cell infiltration and characterizing the functional dimension of immune dysregulation operating in the synovial membrane (Yu et al., 2021).

Data acquisition

Gene expression datasets (GSE55457, GSE12021, and GSE77298) consisting of synovial tissue samples were retrieved from the GEO repository (https://www.ncbi.nlm.nih.gov/), a widely used open-access platform (Edgar et al., 2002; Clough et al., 2024). Corresponding platform-specific annotation files were also saved for further probe-to-gene mapping. An additional independent dataset (GSE89408), also derived from synovial tissue samples, was used as an independent validation cohort to further evaluate the identified hub genes (Clark et al., 2012). A detailed characteristic of all the GEO datasets is provided in Table 1.

Data preprocessing and batch-effect correction

By means of the GEOquery package in R, the GEO Series Matrix files for three training microarray datasets, GSE12021, GSE55457, and GSE77298, were downloaded. Each dataset was independently processed. The probe IDs were converted to gene symbols using the annotation files of the particular platforms. Unassigned genes that do not have valid gene symbols were removed. For duplicate probes mapping to the same gene, the expression values were averaged using the avereps() function of the limma package in R (Yang et al., 2022). Only RA and healthy control (HC) samples were chosen for further analysis. The expression datasets of each particular dataset were preprocessed according to the above procedure. After probe annotation, probes without valid gene symbols were removed, and duplicate probes were summarized using the avereps() function. No additional low-expression filtering was made, as the GEO Series Matrix files were preprocessed and normalized by the authors of the original studies. A study used the combined dataset from different publications and data sources. Prior to differential expression analysis, the common gene expression matrix was corrected for batch effects using the ComBat function from the sva R package, with study as a batch variable (Diboun et al., 2006). Principal component analysis (PCA) was performed before and after batch-effect correction to assess the effect of batch correction on the overall expression structure of the samples. The expression matrix after batch-effect correction is used for all further analysis (Chen et al., 2005).

Identification and visualization of differentially expressed genes

DEGs between RA and control synovial tissue samples were identified using the limma package (Ritch et al., 2015). A design matrix (RA versus control) was built, linear models were fitted to the expression data of all samples via the lmFit() function, and eBayes moderation was applied to estimate gene-wise variances (Smyth et al., 2005). DEGs were selected using the criterion of the log2 fold change and adjusted p -value according to the Benjamini–Hochberg method. The genes with |log2FC|>1 and adjusted p -value < 0.05 were considered as DEGs. Enrichment and network analysis were carried out using DEGs (Li et al., 2025). The volcano plot was created using the ggplot2 (Wickham, 2016) package, where log2 fold change (log 2 FC) was plotted along the x-axis and negative log 10 (adjusted p -value) was plotted along the y-axis. Genes with |log 2 FC|>1 and adjusted p -value < 0.05 were used as the criteria to select DEGs with high differences in gene expression. The heatmap of the top-50 DEGs was created using the pheatmap package (Kolde, 2010).

Functional enrichment analysis

Functional enrichment of DEGs was carried out using the clusterProfiler package in R (Yu, 2024). Gene symbols of the DEGs were converted into Entrez Gene IDs for subsequent enrichment analysis (Inamo, 2020). We identified enriched DEGs in biological process, molecular function, and cellular component. GO enrichment analysis was performed for the biological process, molecular function, and cellular component terms. KEGG pathway enrichment analysis was also conducted to observe enriched DEGs. The Benjamini–Hochberg method was used for multiple testing correction, with adjusted p -value <0.05 considered statistically significant for GO and KEGG enrichment analyses. For GO and KEGG enrichment analyses, gene ratio and the corresponding gene count (Count) were reported for each enriched term or pathway. Gene set enrichment analysis (GSEA) was performed using the fgsea algorithm and the gseKEGG() function from the clusterProfiler package in R, which uses the KEGG database for gene sets. This method was applied to identify the biological characteristics and functions of all the genes in the training set. Genes were ordered based on their logFC scores in a descending order, and 1,000 permutations were run for GSEA (Subramanian et al., 2005; Reimand et al., 2019).

Protein–protein interaction network analysis

To explore potential connections among DEGs, a protein–protein interaction (PPI) network was built. The top-100 DEGs, ranked by absolute logFC, were submitted to the STRING database with a combined interaction score cutoff of ≥0.400 (Szklarczyk et al., 2023), which was viewed and analyzed using Cytoscape. The Cytoscape plugin’s maximum clique centrality (MCC) algorithm was applied to find hub genes (Shannon et al., 2003; Chin et al., 2014).

Machine learning analysis

To refine the candidate hub genes into diagnostic biomarkers, machine learning algorithms were applied (Wu et al., 2024). The training dataset consisted of the combined microarray data that include 41 RA and 26 HC samples. A total of 20 hub genes from the PPI network were selected as the feature set, and the disease status (RA vs. HC) was selected as the response variable.

LASSO (least absolute shrinkage and selection operator) logistic regression was carried out using the glmnet package in R (Friedman et al., 2010). A binomial logistic regression model was fitted. To find the optimum regularization parameter (λ), stratified 10-fold cross-validation was performed. Candidate biomarkers were selected based on non-zero regression coefficients; the final model was created using the lambda.1se criterion (λ = 0.1015).

Support vector machine-recursive feature elimination (SVM-RFE) was carried out using the caret package along with a linear support vector machine classifier (Kuhn, 2008). Stratified five-fold cross-validation repeated 10 times was performed for all feature sets that contained 1–20 hub genes. The feature set with the highest cross-validation accuracy was chosen as the optimal one. The final hub genes were the genes that overlapped with LASSO and SVM-RFE. All feature selection and modeling were only performed with the training dataset (41 RA and 26 HC). The independent RNA-seq dataset GSE89408 (152 RA and 28 HC) was not involved in feature selection and modeling procedures.

External validation

The potential hub genes obtained in this study were subjected to external validation through an independent RNA-seq dataset, namely, GSE89408. The dataset included 152 RA and 28 HC samples. No class-balancing procedure was applied because this dataset was used only for independent validation. The unequal distribution of RA and HC samples may influence classification-derived performance measures, particularly accuracy, sensitivity, specificity, and F1 score; therefore, the validation results were interpreted primarily based on the ROC-AUC values and were considered supportive rather than definitive evidence of diagnostic performance. Raw read counts were normalized using the edgeR package (Chen et al., 2025) by estimating normalization factors and transforming the results into log2 counts per million (CPM). The normalized expression levels of CXCL13 and CXCL10 in RA and HC samples were then analyzed for a consistent expression pattern.

ROC curve analysis

The ROC curve analysis was carried out using the pROC package in R (Robin et al., 2011) to determine the diagnostic potential of the detected hub genes. ROC curve was drawn for the CXCL13 and CXCL10 genes. The area under the ROC curve (AUC) was used to evaluate the diagnostic potential of the models.

Immune cell infiltration analysis

The relative abundance of 22 immune-cell subsets in RA and healthy control synovium was estimated with CIBERSORTx (https://cibersort.stanford.edu/) (Newman et al., 2019) using the LM22 signature matrix and 1,000 permutations. The batch-corrected microarray expression matrix was used as the input, with quantile normalization enabled and batch correction disabled. The analysis was performed in the relative mode using the log2-transformed batch-corrected microarray expression matrix. All 67 samples satisfied the CIBERSORTx deconvolution quality criterion (p < 0.05). Immune cells in the two groups were compared using the Wilcoxon rank sum test, and the resulting p -values were adjusted for multiple testing using the Benjamini–Hochberg method. Next, Spearman correlation analysis was performed between the expression levels of CXCL13 and CXCL10 and the estimated fractions of all 22 immune-cell subsets in the RA samples.

Data integration and batch correction

The three GEO datasets were merged. After annotating probes, eliminating those without valid gene symbols, and summarizing duplicate probes, GSE77298, GSE55457, and GSE12021 retained 22,880, 13,237, and 13,237 unique genes, respectively (Supplementary Table S1), and the flow chart is shown in Supplementary Figure S1. The combined expression matrix included 13,237 common genes for subsequent analyses. To validate the hub genes identified, the RNA-seq dataset GSE89408 served as an independent validation cohort. Figure 1A illustrates that the PCA plot of the combined data revealed a significant batch effect, as the three datasets were separated by study. Post-ComBat normalization, the data demonstrated better sample integration across the studies (B). The distinction between RA and control samples remained evident when the samples were categorized by disease status (RA and healthy control), both prior to (C) and following batch-effect correction (D). The empirical Bayes approach successfully removes batch effect without having much impact on biological variance. Batch-effect correction successfully removed the majority of technical variations, thus leading to greater similarity of samples from different batches while maintaining the biological separation between RA and the control samples. This indicates that ComBat effectively minimized technical variability while maintaining the biological signal associated with the disease. Subsequent analyses utilized the expression data that had been normalized and adjusted for batch effects.

Identification of differentially expressed genes

The differential expression analysis conducted using the limma package led to 440 DEGs in the comparison of RA and healthy control samples with |log2FC| > 1 and adjusted p -value <0.05 after multiple test corrections. Specifically, 243 upregulated genes and 197 downregulated genes were identified in RA synovial tissue. Genes associated with immune response processes were mostly upregulated (for example, chemokines and cytokine receptors), while genes related to metabolism and structural proteins were mostly downregulated. A volcano plot was generated using the ggplot2 package to visualize the distribution of DEGs, as illustrated in Figure 2 below. The heatmap of the top-50 DEGs was generated using the pheatmap package, shown in Figure 3. Gene expression data were scaled, and the genes and samples were hierarchically clustered. As shown in the picture below, the gene expression patterns were clearly distinct between the RA and control samples, with some clustering within each group. Several gene clusters were either upregulated or downregulated in RA samples as compared with that in the control ones, indicating an association of these genes in disease-specific gene expression patterns.

Functional enrichment analysis

The DEGs were significantly enriched in immune-related biological processes, including inflammatory response, T-cell activation, and B-cell receptor signaling pathway. The complete GO enrichment results, including the gene ratio and gene count, are provided in Supplementary Table S2. The proteins predicted to be encoded by DEGs in this category are associated with extracellular signal transduction and immunity. Figure 4. The DEGs were notably enriched in immune-related pathways, including cytokine signaling, cytokine–cytokine receptor interaction, and chemokine-mediated signaling, which are closely associated with RA, as shown by KEGG pathway analysis (Figure 5). The complete KEGG enrichment results, including the p.adj value and gene count, are provided in Supplementary Table S3. GSEA based on the ranked logFC gene list identified several significantly enriched KEGG pathways related to the pathogenesis of RA. The significantly enriched pathways were chosen based on the criteria of P -value <0.05 and FDR (q-value) < 0.25, and the top enriched pathways are shown in Figure 6. The complete GSEA results, including gene-set size, normalized enrichment score (NES), p-value, FDR-adjusted p-value, and ranks, are provided in Supplementary Table S4.

PPI network and hub gene identification

The PPI network was constructed using the top-100 DEGs. The maximum clique centrality (MCC) algorithm implemented in the cytoHubba plugin of Cytoscape was applied to identify hub genes. The top 20 genes ranked by the MCC algorithm were identified as hub genes (Figure 7). The genes are mostly associated with immune cells, chemokine signaling, and inflammatory responses, indicating their crucial roles in RA pathogenesis.

Machine learning approaches

Machine learning methods were utilized for the identification of optimal biomarkers from the set of candidate genes derived from the PPI network. LASSO logistic regression analysis was performed for 20 hub genes with the 10-fold cross-validation method. The cross-validation plot depicting the choice of the penalty parameter (λ) is presented in Figure 8A. According to the lambda.1se criterion (λ = 0.1015), only two genes (CXCL13 and CXCL10) having non-zero regression coefficients were kept by LASSO (Supplementary Table S5). Coefficients of the candidate genes at all possible values of λ are depicted in Figure 8B, and it can be seen from the figure that the coefficients of most of the candidate genes shrink to zero with increasing value of λ. Based on the chosen optimal λ from the cross-validation plot (Figure 8A), the final LASSO model was decided, and genes having non-zero coefficients at λ were selected as the optimal LASSO-selected biomarkers.

SVM-RFE was performed by a linear SVM classifier with five-fold repeated 10-times cross-validation. Different feature subsets comprising 1–20 hub genes were analyzed, and the optimal set of genes was selected on the basis of maximum cross-validation accuracy (Figure 9). CXCL13 and CXCL10 were selected as the optimal biomarkers according to SVM-RFE analysis (Supplementary Table S6).

Identification of final hub genes

The final hub genes were the intersection set of the selected genes from the two machine learning algorithms utilized. CXCL13 and CXCL10 were identified by both the LASSO and SVM-RFE algorithms; they are the potential candidate diagnostic biomarkers of the present study. These biomarkers were then evaluated in an independent RNA-seq dataset GSE89408.

ROC analysis of candidate biomarkers

After ROC curve analysis of the identified biomarkers, CXCL13 produced an AUC value of 0.954, an accuracy value of 89.6%, a sensitivity of 85.4%, a specificity of 96.2%, and an F1 score of 90.9% in the training set. CXCL10 produced an AUC value of 0.925, an accuracy value of 88.1%, a sensitivity of 87.8%, a specificity of 88.5%, and an F1 score of 90.0% in the training set. The ROC curve of the biomarkers in the training set is presented in Figure 10, while the expression profile is presented in Figure 11.

After evaluating the biomarkers using the independent RNA-seq validation dataset (GSE89408), CXCL13 produced an AUC of 0.901, an accuracy of 81.7%, a sensitivity of 80.9%, a specificity of 85.7%, and an F1 score of 88.2%. CXCL10 produced an AUC of 0.886, an accuracy of 78.9%, a sensitivity of 75.7%, a specificity of 96.4%, and an F1 score of 85.8%. The ROC curve of the biomarkers in the validation set is presented in Figure 12, while the expression profile is presented in Figure 13.

The AUC values of the biomarkers in both the training and validation datasets are summarized in Table 2. Although the performance of the biomarkers was lower in the validation set than in the training set, they still retained good diagnostic performance. This shows the potential use of CXCL13 and CXCL10 as diagnostic biomarkers for RA. A detailed summary of the performance metrics, including accuracy, sensitivity, specificity, precision, and F1 score for both datasets, is provided in Supplementary Table S7.

Immune cell infiltration in RA synovium

Figure 14 depicts the immune cell infiltration profile analyzed by CIBERSORTx based on the LM22 signature in the samples of RA patients and healthy individuals. Immune cell infiltration in the synovium of RA patients and healthy individuals differed. Compared to healthy subjects, the synovium of RA patients showed significant differences in plasma cells, T cells follicular helper, gamma delta T cells, monocytes, M1 macrophages, as well as decreased infiltration of resting mast cells. The mean immune cell fractions for the RA and control groups are presented in Supplementary Table S10, while the significantly altered immune cell populations are summarized in Supplementary Table S9. Importantly, the hub genes identified exhibited significant associations with different immune cell infiltrates analyzed by CIBERSORTx. CXCL13 was associated with multiple immune-cell populations, including plasma cells, T-cell subsets, and M1 macrophages, whereas CXCL10 showed significant associations with M1 macrophages, plasma cells, and selected T-cell subsets. These differences in immune cell infiltration between the two groups are illustrated in Figure 15.

Hub genes and immune cells correlation analysis

Correlation analysis between RA-associated hub gene expression and immune cell infiltration was conducted on the RA cohort using estimates of the fractions of 22 immune cell populations and the expression of CXCL10 and CXCL13 (Figure 16). Multiple testing correction using the Benjamini–Hochberg procedure identified significant positive correlations between CXCL13 expression and M1 macrophages (r = 0.75) and between CXCL10 expression and M1 macrophages (r = 0.83), and these along with complete correlation results across all 22 immune cell populations are provided in Supplementary Table S11. The strongest significant correlations are presented in Supplementary Table S8 and summarized in Table 3.

RA is a long-term inflammatory condition that targets the immune system and produces discomfort in the joints in the hands and legs. In the initial stage, the joint lining (synovial membrane) is affected. In addition, despite advances in understanding the molecular mechanisms underlying the disease, identification of reliable early diagnostic biomarkers remains limited, mainly due to disease heterogeneity (Paleolog, 2002). Moreover, we used integrative bioinformatics and machine learning approaches to generate an analytical platform that can be used to identify and validate potential biomarkers for the diagnosis of RA and elucidate their associated biological pathways (Yang et al., 2022).

After combining the gene expression data, 440 genes were significantly identified. The application of batch-effect correction before differential expression analysis improved the integration of the three independent datasets while preserving the biological differences between RA and HC samples, thereby increasing the reliability of the downstream analyses. By combining gene expression data from RA and control samples, the results revealed greater stability and reduced variance, including across experimental settings. Additionally, DEGs in RA and control samples are clearly distinguishable, including gene expression changes in RA. By using functional enrichment analysis methods (GO, KEGG, and GSEA) (Reimand et al., 2019), we also identified that DEGs were enriched in the chemokine signaling pathway, cytokine–cytokine receptor interaction, primary immunodeficiency, and hematopoietic cell lineage. By constructing and analyzing the PPI network (Shannon et al., 2003), we identified 20 genes that were significantly involved in gene interactions. Chemokine activity and immune cell activation were the main functions of these genes. Two different machine learning algorithms (LASSO and SVM-RFE) (Friedman et al., 2010; Kuhn, 2008) were applied to improve the robustness of feature selection and minimize the risk of overfitting. The use of two complementary feature-selection algorithms reduced the likelihood of selecting dataset-specific biomarkers and increased confidence in the final candidate genes identified by both methods. As a result, the hub genes CXCL13 and CXCL10 remained in the network. An independent test dataset was used to validate the expression of these genes in RA. ROC curve analysis of these genes indicated that they are promising biomarkers for RA diagnosis (Robin et al., 2011).

Furthermore, we deployed CIBERSORTx (Newman et al., 2019) to determine immune cell infiltration in RA patients’ synovial tissues. Compared with the controls, RA synovial tissues show increased numbers of plasma cells, M1 macrophages, and γδ T cells but decreased numbers of resting mast cells, monocytes, and NK cells. The correlation between the hub genes and immune cells was verified by correlation analysis (Lu et al., 2025). CXCL13 and CXCL10 had a significant association with activated immune cells, such as macrophages and plasma cells, but inversely correlated with resting immune cells. These results indicate that CXCL13 and CXCL10 are associated with immune cell infiltration patterns. CXCL13 has been reported to be associated with B-cell migration and is also involved in the formation of ectopic lymphoid structures in the synovial tissue. CXCL10 has been reported to be associated with T-cell migration and macrophages during inflammation. Meanwhile, these genes were also closely related to the movement of leukocytes, cytokine signal transduction, and interactions between leukocytes, which were all closely related to the pathogenesis of RA (Choy, 2012; Zlotnik and Yoshie, 2000).

The CXCL13–CXCR5 and CXCL10–CXCR3 signaling axes may represent potential target pathways that warrant further investigation. The CXCL13–CXCR5 pathway is important for B-cell migration (Pan et al., 2022). The CXCL10–CXCR3 pathway is important for the migration of T cells and macrophages (Yue et al., 2024). Inhibiting CXCR3 and CXCR5 receptors on immune cells may also block activation by other ligands within the multi-ligand, multi-receptor families of chemokines. This study is based on the analysis of gene expression using microarray analysis, which may have inherent limitations. We found that the levels of the two chemokines CXCL13 and CXCL10 were high in patients with RA and associated with immune dysfunction. As these two proteins are highly relevant and valuable for the diagnosis of RA, their corresponding receptor axes (CXCR5 and CXCR3) may represent potential hypotheses for future investigation (Zhao et al., 2021). It is worth noting that CXCL13 and CXCL10 have previously been examined in RA, including as baseline serum predictors of response to TNF-inhibitor therapy in a small prospective clinical cohort (Han et al., 2016). The present study differs in both the aim and approach: rather than a treatment-response biomarker measured in peripheral blood, we establish CXCL13 and CXCL10 as diagnostic biomarkers derived from synovial tissue transcriptomics, using an integrative framework that combines batch-corrected multi-cohort gene expression data, dual machine-learning feature selection (LASSO and SVM-RFE), and external validation in an independent RNA-seq dataset. We further link these biomarkers to specific immune cell infiltration patterns via CIBERSORTx deconvolution, providing a mechanistic dimension that prior single-cohort protein studies did not address. Overall, this work highlighted CXCL13 and CXCL10 as possible biomarkers for diagnostic and potential future exploration of RA.

This research explored the possible novel biomarkers and pathways for RA using bioinformatics and machine learning methodologies. DEGs were first identified by integrating various gene expression patterns and reducing batch effects (Wu et al., 2024). DEGs were implicated in immune-related biological processes, including the interferon signaling pathway, inflammatory responses, and cytokine–cytokine interaction, according to functional enrichment analysis, indicating a strong relationship between immunological dysfunction and the development of RA. According to protein–protein interaction network analysis and machine learning (LASSO and SVM-RFE) techniques (Wu et al., 2024), CXCL13 and CXCL10 were hub genes that were upregulated in RA synovial tissues and showed a high area under the ROC curve (AUC).

We identified the biological processes linked to the gene expression profile using immune infiltration analysis, including a pro-inflammatory state due to high concentrations of activated immune cells, including macrophages and plasma cells. The potential involvement of chemokines such as CXCL13 and CXCL10 in RA was assessed using correlation analysis (Han et al., 2016). However, the receptor pathways through CXCR5 and CXCR3 may be studied as potential targets for future research (Zhao et al., 2021).