Deciphering the molecular landscape of rheumatoid arthritis offers new insights into the stratified treatment for the condition

Abstract

Background: For Rheumatoid Arthritis (RA), a long-term chronic illness, it is essential to identify and describe patient subtypes with comparable goal status and molecular biomarkers. This study aims to develop and validate a new subtyping scheme that integrates genome-scale transcriptomic profiles of RA peripheral blood genes, providing a fresh perspective for stratified treatments.

Methods: We utilized independent microarray datasets of RA peripheral blood mononuclear cells (PBMCs). Up-regulated differentially expressed genes (DEGs) were subjected to functional enrichment analysis. Unsupervised cluster analysis was then employed to identify RA peripheral blood gene expression-driven subtypes. We defined three distinct clustering subtypes based on the identified 404 up-regulated DEGs.

Results: Subtype A, named NE-driving, was enriched in pathways related to neutrophil activation and responses to bacteria. Subtype B, termed interferon-driving (IFN-driving), exhibited abundant B cells and showed increased expression of transcripts involved in IFN signaling and defense responses to viruses. In Subtype C, an enrichment of CD8+ T-cells was found, ultimately defining it as CD8+ T-cells-driving. The RA subtyping scheme was validated using the XGBoost machine learning algorithm. We also evaluated the therapeutic outcomes of biological disease-modifying anti-rheumatic drugs.

Conclusions: The findings provide valuable insights for deep stratification, enabling the design of molecular diagnosis and serving as a reference for stratified therapy in RA patients in the future.

Keywords: gene expression profiles; machine learning; rheumatoid arthritis; stratification; unsupervised clustering.

Figure 1

The workflow of data processing procedures in the study. Four microarray datasets containing 1,138 RA patients and three RNA-seq datasets including 268 RA patients were selected as training sets and test sets respectively from the public database. DEGs were filtered after the normalization, and unsupervised clustering was performed with enrichment analysis followed. Then, the XGBoost algorithm was contrived to predict the responses of stratified subtypes to commonly used five treatments. Finally, in search for more novel therapeutic targets of RA patients, drug prediction was carried out by utilizing SMR and the Open Targets platform. DEGs, differentially expressed genes; SMR, Summary data-based Mendelian randomization analysis.

Figure 2

Identification of differentially expressed genes (DEGs) between patients with rheumatoid arthritis (RA) and healthy controls (HCs). (A, B) The heatmap and volcano plot of DEGs in RA patients versus HCs. (C–E) GO enrichment, KEGG and Reactome analyses of 163 up-regulated DEGs.

Figure 3

Identification and gene expression characterization of rheumatoid arthritis (RA) subtypes. (A) The consensus score matrix for RA samples when k = 3. A higher consensus score between the two samples indicated they were more likely to be assigned to the same cluster in different iterations. (B) Consensus clustering for the cumulative distribution function for k = 2–6. (C) Relative changes in the area under the cumulative distribution function curve for k = 2–6. (D–F) Molecular pattern distribution of three subtypes of RA in different biological processes and pathways. The top 20 most significantly enriched biological processes in each subtype of GO BP database and the top 5 most important signaling pathways in the Reactome database. (G) A Venn diagram showing up-regulated DEGs in subtype A, subtype B and subtype C compared with HCs.

Figure 4

Pathway and cell subset-driven characterization in RA subtypes. (A, B) Enrichment scores for pathways and cell subsets for each RA subtype. Box plots for the enrichment scores of pathways and cell subsets for each RA subtype. Wilcoxon test was used to analyze the differences across three subtypes. ns, not significant; *P<0.05; **P<0.01; ***P<0.001; **** P<0.0001. FDR, false discovery rate.

Figure 5

Distribution of gene-driven subtypes and multiple biologics treatments respond to the RA subtypes. (A–C) The variations across the three subtypes concerning autoantibodies, disease activity and four immune cells. The box plots show the disease activity scores of the three subtypes and the enrichment scores of immune cells. Responder: responded to the biologics; non-responder: did not respond to the biologics. (D) Responder/non-responder to infliximab: 29.5%/70.5% in subtype A, 40.9%/59.1% in subtype B and 46.7%/53.3% in subtype C. (E) Responder/non-responder to Tocilizumab: 22.2%/77.8% in subtype A, 35.3%/64.7% in subtype B and 36.4%/63.6% (0/5) in subtype C. (F) Responder/non-responder Rituximab: 59.1%/40.9% in subtype A, 63.6%/36.4% in subtype B and 72.0%/28.0% in subtype C. (G) Responder/non-responder to Abatacept 14.3%/85.7% in subtype A, 21.4%/78.6% in subtype B and 30.0%/70.0% in subtype C. (H) Responder/non-responder MTX: 40.0%/60.0% in subtype A, 62.5%/37.5% in subtype B and 37.5%/62.5% in subtype C. Wilcoxon test was used to analyze the differences across three subtypes. ns, not significant; *P<0.05; **P<0.01; ***P<0.001; **** P<0.0001.

Figure 6

Manhattan plot for the 9 proteins identified in RA. Each point in the plot indicates a single association test between a plasma protein and RA as the -log10 (P) of a z-score test result which is ordered by genomic position on the x-axis and the association strength on the y-axis. The red horizontal line represents the significant threshold for the P value of FDR less than 0.05 under Bonferroni correction.