Deciphering Immunometabolic Landscape in Rheumatoid Arthritis: Integrative Multiomics, Explainable Machine Learning and Experimental Validation

Introduction

Rheumatoid arthritis (RA) is a chronic, systemic autoimmune disease characterized by persistent polyarticular inflammation, synovial hyperplasia, pain, joint swelling and stiffness.¹ The global prevalence of RA continues to rise as the population ageing, with the number of people living with the disease has more than doubled in the last three decades and was predicted to reach 317 million worldwide by 2050.² As the disease progresses, RA can lead to progressive joint destruction and deformity, resulting in long-term disability, chronic pain and even premature death. Notably, even though several therapeutic strategies, such as biologics and small molecule medications, have proven effective for RA, about 40% of patients do not achieve clinical remission.³ Additionally, RA could increase the risk of infection, depression, malignancy, cardiovascular and respiratory disease.^4,5 The impending increase underscores the urgent need to amplify research endeavors, further exploration of pathogenesis, and improved treatment strategies to mitigate the impact of this prevalent autoimmune disease.

The complex immunological mechanisms of RA have been a focal point of research, with the aberrant activation and infiltration of macrophages within the synovium, a hallmark feature that contributes to the relentless progression of joint damage. The clinical assessment of RA is often guided by the peripheral blood leukocyte count, a critical marker that reflects the disease’s activity and severity. The Th17 subset of CD4+ T cells and their signature cytokine, interleukin-17 (IL-17), have emerged as pivotal mediators in the inflammatory cascade of RA, exacerbating synovial inflammation and hastening cartilage degradation.⁶ Furthermore, recent advancements in immunological research have highlighted the critical involvement of dendritic cells and T cells in the progression of RA.⁷ Dendritic cells (DC), as antigen-presenting cells, are instrumental in shaping the adaptive immune response, while T cells orchestrate a series of downstream inflammatory events. Their interplay is thought to be crucial in the initiation and perpetuation of RA, offering novel targets for therapeutic intervention.

The inappropriate growth of immune and stromal cells, leading to a metabolic disarray, forms the cornerstone of the complex interplay between immune function and inflammation in the disease.⁸ The distinctive metabolic signatures in T cells and macrophages, the metabolic status of lymphocytes, and the metabolic dysregulations in immune and synovial cells are all intricately woven into the tapestry of RA development.^9,10 Meanwhile, adiponectin and leptin secreted by adipose tissue correlate with disease duration and activity, with lipocalin slowing symptoms while leptin promoting inflammation.¹¹ Moreover, adenosine deaminase (ADA) has been implicated in the immunological reshaping of RA, modulating immune responses through the degradation of adenosine, thereby contributing to the progressive joint destruction.¹² The relationship between immunity and metabolism, termed immunometabolism, is bidirectional and complex. Research has established that immunometabolism is not only a byproduct of the inflammatory cascade but also contributes to the initiation and perpetuation of inflammation.^13,14 However, the precise mechanisms underlying this involvement remain elusive. Therefore, exploring the complex immune metabolic landscape in RA is critical for identifying new and effective therapeutic strategies that have the potential to halt or reverse the debilitating effects of the disease.

In this study, we adopt a comprehensive approach to elucidate the complex immunometabolic landscape of RA. Initially, we employ Mendelian randomization (MR) and mediation analysis to investigate the underlying causal relationships between immunity, metabolism and RA. This is followed by single-cell RNA sequencing (scRNA-seq) to explore the role of specific metabolites in immune cells, with the aim of identifying unique gene sets that could serve as novel therapeutic targets. We then apply three advanced machine learning (ML) algorithms to identify biomarkers with diagnostic and therapeutic significance by large bulk sequence data. Utilizing Shapley Additive Explanations (SHAP), an explainable AI tool, we enhance the interpretability of algorithms while preserving the prediction accuracy of complex machine learning models, offering a clear understanding of the model predictions.¹⁵ Finally, we validate these promising biomarkers using clinical samples to ensure their real-world relevance. This integrative approach will facilitate a deeper understanding of the molecular mechanisms linking metabolism and RA, paving the way for the development of precision medicine strategies tailored to the immunometabolic profiles of individual patients.

Material and Methods

Study Design and Data Sources

The study protocol is presented in Figure 1. First, we obtained published GWAS summary data that included traits such as plasma metabolites, peripheral immune cells, and RA. Based on those data, we performed a bidirectional two-sample MR to explore potential causal relationship. Second, two-step MR analysis was used to identify the mediation effect of plasma metabolites on the relationship among metabolites, immune cells, and RA. The MR study was conducted in accordance with the STROBE-MR guidelines (Table S1).¹⁶ Third, we tried to elucidate the variances in metabolite-associated biological pathways within immune cells at the scRNA-seq level. Additionally, we investigated the specific contributions of immune cells exhibiting distinct metabolic profiles throughout the disease progression. Fourth, we have employed three advanced boosting machine learning algorithms to identify pivotal differentially expressed genes (DEGs) at the transcriptomic level. Finally, we procured a series of clinical specimens and corroborated the expression patterns of the core mRNA transcripts through quantitative real-time polymerase chain reaction (qRT-PCR) analysis. Detailed information about the data sources has been summarized in the Table S2 and all the data were publicly available.

Figure 1 Overview study of the study design. (A) Bidirectional MR was implemented to investigate the causal relationships among plasma metabolome, immune cells and RA. (B) Two-step MR was utilized to estimate the mediation effect of immune cells in the relationship between plasma metabolome and RA. (C) scRNA-seq was employed to distinguish the disparities in DC derived from CT and RA. To identify these variations, we performed a comprehensive analysis using bulk RNA-seq data. Subsequently, machine learning techniques were applied to pinpoint key genetic determinants. Abbreviations: MR, Mendelian randomization; RA, rheumatoid arthritis; DC, dendritic cell; CT, control.

Mendelian Randomization Analysis

The genetic variants used as instrumental variables (IVs) should meet three basic assumptions: 1) IVs exhibit a correlation with the exposure; 2) IVs are not associated with confounders; 3) IVs affect the outcome only via the exposure.¹⁷ The single nucleotide polymorphisms (SNPs) were selected based on their significant association with plasma metabolites (P < 1×10⁻⁵) and immune cells (P < 5×10⁻⁶).^18,19 To avoid the confounding effects of linkage disequilibrium (LD), we applied an exclusion criterion with an R² threshold of 0.001 and a maximum distance threshold of 10,000kb.²⁰ If SNPs were associated with the same exposure, we prioritized the SNP with the lowest P-value for further analysis. Additionally, SNPs with an F-statistic below 10 were excluded from our analysis, as they are considered weak instruments that could introduce bias into the results.²¹ We considered an IV to exhibit horizontal pleiotropy if the MR-Egger regression intercept was not equal to 0 and had a statistically significant P-value less than 0.05. In the reverse MR analysis, a P-value greater than 0.05 indicated no evidence of reverse causality between the exposure factors and the outcome variables.

Mediation Analysis

We applied a two-step MR approach to calculate mediating effects using: Mediated Effect = β(A) × β(B), where β(A) is the effect of exposure on the mediator, and β(B) is the effect of the mediator on the outcome.²² The total effect of metabolites on RA was derived from previous two-sample MR, with the direct effect calculated as Total Effect - Mediated Effect. The proportion mediated was calculated as (Mediated Effect / Total Effect) × 100%, with confidence intervals (CIs) obtained using the delta method. Criteria for mediating roles of immune cells included causal relationships between exposure, mediator, and outcome, with a significant mediated effect. This rigorous methodology ensures robust findings and insights into the immunometabolic pathways linking plasma metabolites to RA development.

Single Cell and Bulk Transcriptome Analysis

In this study, we used the “Seurat” package (version 4.3) to analyze scRNA-seq data from healthy controls (CTs) and RA patients. Quality control excluded cells with fewer than 500 detected genes, over 15% mitochondrial gene content, or over 1% hemoglobin gene content, resulting in 16,482 cells for analysis. We performed principal component analysis (PCA) on the top 2,000 variable genes, addressed batch effects with the harmony algorithm, and used uniform manifold approximation and projection (UMAP) for visualization. Cell types were annotated based on specific marker expression. We calculated gene set scores with the ‘AddModuleScore’ function to reclassify subpopulations, used gene set enrichment analysis (GSEA) to identify biological process differences, and estimated transcription factor activity with the ‘dorothea’ package. We employed the ‘CellCall’ R package to analyze cell-cell communication, providing insights into interactions between cellular subpopulations and the immunological landscape in health and RA.²³

We obtained bulk transcriptome datasets GSE93272 and GSE68689 from the GEO database for training and validation. Using the “Affy” package, we processed the raw data and applied the “ComBat” function from the “sva” package to correct batch effects. DEGs analysis was conducted with the “limma” package, followed by gene set variation analysis (GSVA) to assess gene set enrichment and activity. We also used the IOBR package for immune infiltration analysis, quantifying 10 distinct immune cell types in the tumor microenvironment.

Explainable Machine Learning Analysis

We employed the least absolute shrinkage and selection separator (LASSO) in the training set for dimensionality reduction and DEGs selection. These regularization methods were used to avoid overfitting and improve the interpretability and prediction accuracy of the model. Then, nine machine learning algorithms, including logistic regression (LR), random forest (RF), decision tree (DT), support vector machine (SVM), k-nearest Neighbor (k-NN), gradient boost machine (GBM), categorical boosting (CatBoost), natural gradient boosting (NGBoost), and extreme gradient boosting (XGBoost) were utilized to identify RA-related immune genes. The diagnostic performance of our model was evaluated using receiver operating characteristic (ROC) curves and the area under the curve (AUC). The SHAP (version 0.45.0) package to calculate the SHAP value for each of the features in the best model we selected, which represents the contribution of individual genes to the predictive model.

Peripheral Blood Mononuclear Lymphocyte Collection and Quantitative Real-Time PCR (qRT- PCR)

Our study complies with the Declaration of Helsinki. Our study involved the collection of peripheral blood samples from patients diagnosed with RA and CTs. Prior to sample collection, we obtained informed consent from all participants. The informed consent from the Institutional Review Board of the First Affiliated Hospital of Jinan University (Ethics Approval Number K*-2024-**2) was presented in clear, non-technical language to facilitate understanding, and participants were given ample time to consider their participation and to ask questions. Blood was collected in a manner that minimized discomfort and risk to the participants. Blood was drawn into EDTA tubes and processed within 4 hours to isolate peripheral blood mononuclear cells (PBMCs) using Ficoll-Paque Plus density gradient centrifugation. The mononuclear layer was washed with PBS, and cell viability was confirmed to be over 95% using trypan blue exclusion. PBMCs were either used immediately or cryopreserved in liquid nitrogen for future analysis.

The described process involves extracting total RNA from PBMCs using TRIzol reagent, assessing its integrity and purity, and then reverse transcribing it into cDNA using RT Master Mix (Transgen, Beijing, China). The cDNA is subjected to real-time PCR using SYBR Green qRT‑PCR Master Mix (Transgen, Beijing, China) for accurate gene expression quantification. Amplification and data acquisition are performed under optimized conditions, with β-actin serving as an internal control. Primer sequences are detailed in Table S3.

Drug Repositioning and Molecular Docking

We employed drug repositioning by inputting the expression levels of hub genes into the Connectivity Map (CMAP) database, a widely recognized tool for screening potential therapeutic compounds. Based on the connectivity scores, we selected the top 30 compounds as candidates. From this selection, we focused on the six highest-ranking small molecules for molecular docking with the identified hub genes. The chemical structures of these small molecules were obtained from the PubChem database to ensure accurate representation. We sourced the protein structures of the hub genes from the Protein Data Bank (PDB) and AlphaFold databases. Molecular docking was performed using AutoDock Vina, where lower binding energy indicate higher affinity between the ligands and target proteins.²⁴ The docking result were visualized using PyMOL, enabling a detailed analysis of the interactions between the small molecules and hub genes.

Phenome-Wide Association Analysis

To further evaluate the potential side effects associated with drug targets, we conducted a phenome-wide association study (PheWAS) using the AstraZeneca PheWAS Portal (https://azphewas.com/). This analysis included approximately 15,500 binary and 1,500 continuous phenotypes from about 450,000 exome-sequenced participants in the UK Biobank. Detailed methodology is available in the original publication.²⁵ We applied multiple corrections and set a threshold of p < 1 × 10⁻⁵, the default value in the AstraZeneca PheWAS Portal, to reduce false positives. We applied multiple corrections and considered corrected phenotypes with p < 1×10⁻⁵ as statistically significant based on previous reports.²⁵

Statistical Analysis

The statistical analysis was executed using Python (version 3.11.5) and R (version 4.3.0) software. Data were expressed as mean along with a 95% CI. Wilcoxon or Student’s t-test was employed to compare differences between two groups, while one-way analysis of variance (ANOVA) was used for comparisons among three or more groups. To determine the correlation between variables, Pearson or Spearman correlation test was used. The statistical p-value was calculated using a two-tailed method and reported with three decimal places, and false discovery rate (FDR) was used to adjust the p-value to minimize false positives.

Results

Causality Between Plasma Metabolites, Immune Cells and RA

To explore the causal impact of plasma metabolomes on RA and immune cells on RA, two-sample MR analyses were performed. The threshold for statistical significance was set as an FDR below 0.05, and the quantity of SNPs for each metric was considered to augment the robustness of the genetic IVs. All these genetic IVs met the requirements of LD-independent (r2 < 0.1) and achieved a genome-wide significance level (p < 1×10⁻⁵).

In our examination of plasma metabolomes and immune cells as exposure variables, we employed a multitude of SNPs as instrumental variables to strengthen our analysis. The result indicated a statistically significant association between six plasma metabolomes, five immune cells and the occurrence of RA after eliminating horizontal pleiotropy of variables by MR-PRESSO (Table S4).

Three plasma metabolomes and immune cells were found to play a protective role in RA pathogenesis, including 18 SNPs for Glutamate to glutamine ratio (OR = 0.815, 95% CI = 0.771–0.859, P = 8.73×10⁻⁵, PFDR = 0.024), 33 SNPs for metabolonic lactone sulfate (MSL) levels (OR = 0.931, 95% CI = 0.916–0.946, P = 6.29×10⁻⁵, PFDR = 0.022), 26 SNPs for 1-oleoyl-2-linoleoyl-GPE (PE) levels (OR = 0.929, 95% CI = 0.912–0.945, P = 1.42×10⁻⁴, PFDR = 0.033), 6 SNPs for CD14+ CD16+ monocyte AC (OR = 0.845, 95% CI = 0.809–0.881, P = 7.55×10⁻⁵, PFDR = 0.014), 1 SNP for HEVM on CM CD4+ (OR = 0.816, 95% CI = 0.772–0.859, P = 6.95×10⁻⁵, PFDR = 0.014), and 9 SNPs for HVEM on T cell (OR = 0.881, 95% CI = 0.855–0.907, P = 3.90×10⁻⁵, PFDR = 0.014).

Three plasma metabolomes and two immune cells were found to be risk factors for the pathogenesis of RA, including 16 SNPs in 4-Acetamidobutanoate levels (OR = 1.187, 95% CI = 1.151–1.223, P = 4.53×10⁻⁵, PFDR = 0.021), 15 SNPs in 5-acetylamino-6-formylamino-3-methyluracil (AFMU) levels (OR = 1.093, 95% CI = 1.079–1.107, P = 9.53×10⁻⁸, PFDR = 1.33×10⁻⁴), 26 SNPs for 1-stearoyl-2-docosahexaenoyl-gpc (GPC) (OR = 1.172, 95% CI = 1.145–1.199, P =5.52×10⁻⁷, PFDR = 3.86×10⁻⁴), 5 SNPs for CD28+ CD45RA- CD8br %CD8br (OR = 1.215, 95% CI = 1.175–1.256, P = 3.97×10⁻⁵, PFDR = 0.014), and 9 SNPs for DC absolute count (AC) (OR = 1.139, 95% CI = 1.110–1.168, P = 1.55×10⁻⁴, PFDR = 0.023) (Figure 2). Conversely, when considering RA as the exposure variable, no immune cell nor plasma metabolome achieved statistical prominence in the conventional purview. The causal effects of all plasma metabolomes and immune cells on RA are listed in Table S5. The SNPs used for MR analyses are listed in Table S6.

Figure 2 MR analyses exploring the relationship between PE, DC and RA. (A) Forest plot of bidirectional MR between PE, DC and RA. (B) The scatter plot of the association between PE levels and DC AC. (C) Forest plot show the relationship between PE levels and DC AC. (D) Leave-one-out plot to visualize causal effect of PE levels on DC AC when leaving one SNP out. (E) Funnel plot based on IVW and MR-EGGER regression to assess potential causal effect of the PE levels on DC AC. Abbreviations: RA, rheumatoid arthritis; AC, absolute count; AFMU, 5-acetylamino-6-formylamino-3-methyluracil; CI, confidence interval; DC, dendritic cell; FDR, false discovery rate; GPC, 1-stearoyl-2-docosahexaenoyl-gpc (18:0/22:6); IVW, inverse variance weighted; MSL, metabolomic lactone sulfate; OR, odds ratio; PE, 1-oleoyl-2-linoleoyl-GPE (18:1/18:2); SNP, single nucleotide polymorphism.

Mediation Analysis

To explore the potential mechanisms of RA occurrence and development, we conducted a mediation analysis to identify the causal pathway from plasma metabolites to RA mediated by immune cells. This analysis focused on plasma metabolites and immune cells previously identified by two-sample MR as being associated with RA. Considering that both plasma PE levels and DC are causally related to RA, and that previous literature suggests that plasma PE levels is an intermediate product of autophagy,^26,27 we further explored the relationship between PE and DC by mediation analysis, using plasma PE levels as the target of our study.

The mediation analysis reveals that plasma PE levels exhibit significant protective direct effects (β = −0.102, 95% CI: [−0.184, −0.021]) on DC AC. DC AC also presenting a significant negative mediation effects (β = −0.013, 95% CI: [−0.026, −0.001], P = 0.035) on RA with 17.999% (95% CI: 8.224%, 27.775%) proportion. These results showed the consistent direction of the total and indirect effects, and that the leave-one-out analysis supported the reliable causal relationship in the two-sample MR study of exposure to outcome, exposure to mediator (Table 1; Figure S1; Table S7).

Table 1 The Mediation Effect of PE on RA via DC

scRNA Analysis of PBMCs from CT and RA

We identified seven distinct cellular clusters from PBMC of CT and RA, according to the marker gene expression (Figure 3A and B). Among these, CD4 T cells were found to be the most prevalent and a significant increase in the proportion of DC was observed in the RA compared to the CT (Figure 3C). Subsequently, we scored biological processes associated with PE and observed that autophagy-related scores were markedly elevated in both RA and DC subsets (Figure 3D). This prompted further investigation into the differences between DCs with high and low autophagy scores (Table S8). In high autophagy score DCs, processes such as the regulation of actin filament organization, leukocyte-mediated cytotoxicity, and regulation of actin filament bundle assembly were enhanced, while processes like protein N-linked glycosylation via asparagine, ubiquitin-dependent ERAD pathway, and ERAD pathway were attenuated (Figure 3E, Table S9). Within the DCs exhibiting high autophagy scores, the activity of transcription factors such as HNF1A, LHX2, and POU5F1 was decreased, while the activity of AR, SP1, and USF1 was increased (Figure 3F, Table S10). These findings suggest a complex interplay between autophagy levels, cellular function, and transcriptional regulation in RA. Subsequently, we investigated the functional differences in cellular interactions between the identified clusters. Our observations indicate that DCs with high autophagy scores may contribute to the pathogenesis of RA through processes such as necroptosis and the activation of the Jak-STAT signaling pathway (Figures 3G and 3H). These pathways are known to play crucial roles in immune cell regulation and could provide novel targets for therapeutic intervention in RA.^28,29

Figure 3 scRNA-seq analysis of peripheral blood of CT and RA. (A) Heatmap displaying the marker genes of cell subpopulations. The redder the color, the higher the gene expression. (B) UMAP plots visualizing cell subpopulations. Different subpopulations of cells are composed of individual points. (C) Bar graphs depicting the proportions of cell subpopulations in CT and RA. The numbers and scale graphs represent the percentage of total cells accounted for. (D) Dot plots depicting PE-related biological processes across various cell types. Dot size represents the cell percentage expressing the gene set, with larger dots for higher percentages. Dot shade reflects the mean expression level, with darker shades for higher expression. (E) Bar plot illustrating biological process driven by differentially expressed genes between DC cell groups with varying autophagy scores. It separately shows the top 10 upregulated and downregulated terms, with corrected p-values shown in red. (F) Heat map showing DEGs encoding transcriptional factors. (G) Bubble plot showing the ligand–receptor relationships between DC and other cell types. (H) Chord diagram showing the interaction between DC and other cell types. Abbreviations: CT, control; RA, rheumatoid arthritis; DC, dendritic cell; PE, 1-oleoyl-2-linoleoyl-GPE (18:1/18:2); DEGs, differentially expressed genes.

Bulk Transcriptome Analysis of PBMCs from CT and RA

We obtained transcriptome samples from the GEO database, using GSE93272 as the training cohort and GSE68689 as the validation set, with batch effects corrected (Figure S2). Figure 4A shows DEGs profiles between CT and RA groups in PBMC. GSVA indicated increased autophagy-related scores in DC in RA (Figure 4B). Using lasso regression, we identified 22 key autophagy-associated genes differentiating CT from RA (Figure 4C). Immune infiltration analysis revealed significant associations between these genes and various immune cells, with 14 notably related to DC (Figure 4D). These findings underscore the intricate interplay between autophagy in DC and the immunopathogenesis of RA.

Figure 4 Bulk RNA-seq analysis conducted on peripheral blood mononuclear cells. (A) Heat map representing the DEGs. (B) Bar plot showcasing the autophagy scores, computed via GSVA, contrasting CT and RA. (C) Lasso regression graphs of DEGs employed for genes selection. (D) Heat map of scores displaying the correlation scores between associated genes and immune cell infiltration. Abbreviations: DEGs, differentially expressed genes; GSVA, gene set variation analysis; CT, control; RA, rheumatoid arthritis. *: P < 0.05, **: P < 0.01, ***: P < 0.001.

Identifying Hub Genes Through Explainable Machine Learning and qRT‑PCR Experimentation

We applied 22 key autophagy-associated genes to nine machine learning algorithms to identify key RA-related gene. The validation set (GSE68689) showed AUCs ranging from 0.84 to 0.97. Considering both diagnostic efficacy and learning curve, the K-NN model performed the best with an AUC of 0.94, demonstrating strong predictive ability (Table S11, Figure 5A and Figure S3).

Figure 5 Machine learning and experimental validation in screening key genes for RA. (A) ROC curves representing the performance of nine machine learning algorithms. (B) SHAP bar plot illustrating the contribution of individual genes to the predictive model. (C) Bee swarm plot demonstrating the significance of hub genes in the k-NN model, with each point denoting a sample. (D) The PPI network of hub genes. (E) Bar charts depicting the mRNA expression levels of hub genes. ***P < 0.001, ****P < 0.0001. Abbreviations: RA, rheumatoid arthritis; CatBoost, categorical boosting; NGBoost, natural gradient boosting; SVM, support vector machine; k-NN, k-nearest neighbor; ROC, receiver operating characteristic; AUC, area under curve; SHAP, SHapley additive explanation; XGBoost, extreme gradient boosting; CT, health control.

To interpret the relationship between these genes and RA, we used SHAP analysis. The most relevant gene was PFN1 (SHAP value = +0.06), followed by SRP14 (+0.04), S100A11 (+0.04), SLC25A5 (+0.02), SAP18 (+0.02), and RTN4 (+0.02). Increased expression of PFN1, SRP14, and S100A11 is associated with RA, while decreased expression of SLC25A5, RTN4, and SAP18 is also linked to the disease. Figures 5B and 5C illustrate the contributions of the top nine and twenty genes to the k-NN model, respectively.

We constructed a protein-protein interaction (PPI) network using the STRING database to explore molecular interactions of DC autophagy-associated genes. PFN1 emerged as a central hub, interacting with other key genes (Figure 5D), suggesting its pivotal role in RA pathophysiology. Furthermore, we validated these findings with qRT-PCR, showing significant upregulation of SRP14, S100A11, PFN1, and SAP18 in PBMCs from RA compared to CT (Figure 5E).

Discovery of Potential Drugs

Through CMAP analysis, we identified 30 potential small molecule drugs targeting the immunometabolic landscape of RA, including vismodegib, zonisamide, mirtazapine, clenbuterol, batimastat, and dapivirine (Table S12). To further validate the binding capacity of these small molecules with hub genes, we employed molecular docking techniques to evaluate the stability of their interactions based on calculated free binding energies. Previous studies suggest that a binding energy below −5 kcal/mol indicates good binding, while values below −7 kcal/mol signify strong affinity.³⁰ The molecular docking results demonstrated that six small molecules exhibited moderate binding potential with four hub genes (Figure 6A). We then visualized the interaction of vismodegib, which showed the strongest binding energy with the hub targets. The three-dimensional interaction model revealed that vismodegib can form stable hydrogen bonds with multiple targets (Figure 6B).

Figure 6 Molecular docking results of potential drugs and hub genes. (A) Molecular docking results heat map of potential drugs and hub genes. (B) Molecular docking mode diagram. (C) Binary traits PheWAS association with PFN1. (D) Continuous traits PheWAS association with PFN1. Abbreviations: PheWAS: phenome-wide association study.

To explore the pleiotropic effects of these genes on various human traits, we conducted a PheWAS at the gene level, which could identify associations between gene-determined protein expressions and specific diseases or traits. As illustrated in Figure 6C and D and Figure S4, only the PFN1 gene showed significant correlations with other traits (P < 0.05). The absence of significant associations for the other three drug targets suggests minimal potential side effects from drugs targeting these genes, reinforcing the reliability of our findings. In contrast, PFN1 exhibited a significant inverse correlation with laboratory test results, indicating that RA medications targeting this gene could lead to abnormalities in related tests. This highlights the importance of careful monitoring of health parameters when using such medications in clinical practice.

Discussion

RA is a systemic autoimmune disease characterized by synovitis of joints with elevated inflammatory markers and autoantibodies, including IL-6, RF and other specific antibodies.³¹ Innate and adaptive immune responses play an indispensable role in the pathogenesis of RA. Increasing number of studies have shown that the complex interaction and activation of infiltrating immune cells are key factors in the formation of synovitis and persistent joint damage. Both ACR³² and EULAR³³ recommend a strategy targeting an outcome of remission or low disease activity, Methotrexate is the first-line treatment for RA, and disease–modifying antirheumatic drugs (DMARDs) are recommended to be started as a second-line treatment for RA. Additionally, Sharma et al demonstrated the therapeutic potential of natural products for RA by targeting inflammatory pathways associated with the disease.³⁴ However, only 40% to 60% of patients achieving low disease activity and up to 20% reaching remission after DMARDs administrated,³⁵ even patients in sustained drug-induced remission can experience disease flares, which both risk incremental joint damage and limit physical function and quality of life.³⁶ Effective biomarkers, especially for the early stage, have not been established due to the significant heterogeneity of RA. Early diagnosis and treatment of RA could effectively prevent the disease progression of 90% patients.³⁷ Over 20 years ago, the concept of preclinical autoimmunity was introduced, highlighting the presence of autoantibodies well before disease symptoms appear. This discovery fundamentally changed our understanding of autoimmune diseases, including RA,³⁸ by establishing a clear separation between disease onset and clinical manifestations. RA is heavily influenced by immune cells. Hemmatzadeh et al³⁹ discussed the roles of T and B cells in its pathogenesis. Baker et al³⁶ identified specific cell subpopulations (CD45RO⁺PD1^hi CD4⁺ and CD8⁺ T cells, and CD27⁺CD86⁺CD21⁻ B cells) that predict arthritis flares during drug withdrawal. Sharma et al⁴⁰ found that the hydroethanolic extract of Nyctanthes arbor-tristis may increase B- and T-cell populations. Our research found a causal relationship between RA and five immune cell types, focusing on DCs due to their crucial role in RA pathogenesis. Using single-cell RNA analysis, we observed a significant increase in DCs in RA compared to controls, consistent with previous studies. Canavan et al⁴¹ found that conventional DCs are more prevalent in RA synovial fluid than in peripheral blood, with increased expression of antigen-presenting and co-stimulatory molecules. These cells, when cocultured with T cells, induce proliferation and secretion of IFN-γ, IL-4, and IL-17. Consequently, it is imperative to identify specific plasma metabolites and analyze the pattern of immune cell infiltration associated with RA to improve prognosis. The aim of this study was to identify plasma metabolites specific for RA and to study their relationship with immune cell infiltration and their impact on RA.

To explore the immunometabolic landscape of RA, we conducted a comprehensive analysis involving 1,400 plasma metabolites, 731 immune cells, 91 inflammatory cytokines, and RA itself. We discovered causal relationships between five immune cell types, six plasma metabolites and RA. We are particularly interested in plasma PE levels and DCs, as they play key roles in the pathogenesis of RA. Notably, both are closely linked to cellular autophagy, which significantly influences RA development. Our MR analysis highlighted that higher plasma PE levels are protective against RA, whereas DCs increase the risk. We further examined the interaction between plasma PE levels and DCs through mediation analysis, revealing that elevated PE levels mitigate the risk associated with DCs. Single-cell analysis showed a significant upregulation of autophagy in PE-related processes in RA patients. This prompted us to investigate the link between DCs and autophagy. Using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, we identified DEGs among DC groups with varying autophagy scores. This connection opens new avenues for exploring targeted therapeutic strategies.

Metabolism plays a crucial role in the activation of DCs. Typically, anabolism enhances immunogenicity, while catabolism promotes tolerogenicity.⁴² For the first time, we utilized large-scale GWAS summary datasets of plasma metabolites and RA to conduct a genetic correlation scan. This innovative approach explores the relationship between plasma metabolism and RA, potentially offering new insights into the genetic mechanisms underlying the disease. Human blood metabolism is shaped by genetic factors, complex regulatory mechanisms, and non-genetic elements.⁴³ Evaluating metabolites provides a promising method for identifying the pathophysiological components of complex diseases. As end products of cellular processes, metabolites illuminate the connections between genotypes and phenotypes.^44,45 However, comprehensive studies on the genetic and causal relationships between RA and blood metabolites remain limited. In our study, sc-RNA seq founded that autophagy-related scores were markedly elevated in both RA and DC subsets. Autophagy, itself is a catabolic process. Many inflammatory processes rely on energy supply and metabolic conditions as key factors in tissue damage, contributing to the pathogenesis of RA. It has been pointed out that abnormalities in autophagy can significantly influence cell proliferation and cell loss, further impacting disease progression.⁴⁶

Previous research indicates that 60% of the risk for developing RA is genetic.⁴⁷ Early diagnosis and treatment can prevent or significantly slow joint damage in up to 90% of patients, averting irreversible disability.⁴⁸ Inspired by studies like Ao et al⁴⁹ who identified 679 DEGs enriched in immune response pathways, and Wang et al⁵⁰ who highlighted CD8+ T-cell related genes as diagnostic markers, we conducted our research. Using bulk transcriptome data, we identified key autophagy-related pathogenic genes in RA, with 14 notably related to DCs. We identified a regulatory network where SRP14, S100A11, PFN1, and SAP18 were significantly up-regulated in RA patients, as verified by qRT-PCR results. To our knowledge, few studies have explored the relationship between these autophagy-related genes and RA. Additionally, we developed an explainable artificial intelligence prediction model using machine learning to identify RA onset. SHAP analysis revealed SRP14 as the most critical pathogenic gene influencing RA onset and progression, suggesting it as a core target for designing therapeutic regimens. These key genes enhance predictive accuracy, marking one of the few studies to apply machine learning in this context. This research lays the groundwork for advancing diagnostic tools and therapeutic strategies.

To identify potential drugs for targeting the immunometabolic landscape of RA, we employed CMAP and molecular docking techniques based on the expression levels of hub genes. This approach allowed us to screen 30 promising compounds, with vismodegib showing significant therapeutic potential. Vismodegib is a small molecule inhibitor that targets the Hedgehog signaling pathway, which is essential for embryonic development and tissue repair. However, activated Hedgehog signaling can also lead to abnormal proliferation and activation of fibroblast-like synoviocytes, worsening synovial inflammation and joint destruction in RA.⁵¹ Inhibiting Smoothened, a key component of this pathway, has proven effective in treating RA.⁵² While vismodegib is primarily used for cancer therapy, it presents a strong candidate for RA treatment. Further studies are needed to validate its therapeutic effects and establish the optimal treatment regimen.

This study presents several advantages in constructing the immunometabolic landscape of RA using multi-omics data. First, we introduced mediation analysis into this framework, enabling us to explore causal relationships for the first time. Second, by integrating multi-omics data and utilizing various machine learning algorithms, we developed a robust model that successfully identified hub genes, which were validated with clinical samples. Additionally, we identified potential small molecule drugs and assessed possible side effects of relevant drug targets through drug repositioning, molecular docking, and PheWAS analysis.

However, this study has limitations. First, our analyses focused on plasma metabolite and immune cell data from European populations, highlighting the need for broader ethnic representation. To address this, we included data from diverse ethnic groups for single-cell and bulk RNA sequencing and used Chinese samples for peripheral blood qRT-PCR validation. Although this helps, the small sample size for qRT-PCR validation indicates the need for larger cohort and longitudinal studies, which are challenging due to ethical considerations in clinical sample collection. We are actively working on expanding the sample size and designing further experimental studies to explore these mechanisms in more detail.

Conclusion

In summary, our study clarifies the causal relationships between plasma metabolites, immune cells, and RA, identifying new therapeutic targets and potential drugs. These findings deepen our understanding of RA’s mechanisms. Future studies focusing on autophagy-related pathogenic genes and Vismodegib are anticipated to offer promising avenues for enhanced prevention, diagnosis, and treatment, with significant clinical applications.