Rituximab versus tocilizumab in rheumatoid arthritis: synovial biopsy-based biomarker analysis of the phase 4 R4RA randomized trial

Patients with rheumatoid arthritis (RA) receive highly targeted biologic therapies without previous knowledge of target expression levels in the diseased tissue. Approximately 40% of patients do not respond to individual biologic therapies and 5–20% are refractory to all. In a biopsy-based, precision-medicine, randomized clinical trial in RA (R4RA; n = 164), patients with low/absent synovial B cell molecular signature had a lower response to rituximab (anti-CD20 monoclonal antibody) compared with that to tocilizumab (anti-IL6R monoclonal antibody) although the exact mechanisms of response/nonresponse remain to be established. Here, in-depth histological/molecular analyses of R4RA synovial biopsies identify humoral immune response gene signatures associated with response to rituximab and tocilizumab, and a stromal/fibroblast signature in patients refractory to all medications. Post-treatment changes in synovial gene expression and cell infiltration highlighted divergent effects of rituximab and tocilizumab relating to differing response/nonresponse mechanisms. Using ten-by-tenfold nested cross-validation, we developed machine learning algorithms predictive of response to rituximab (area under the curve (AUC) = 0.74), tocilizumab (AUC = 0.68) and, notably, multidrug resistance (AUC = 0.69). This study supports the notion that disease endotypes, driven by diverse molecular pathology pathways in the diseased tissue, determine diverse clinical and treatment–response phenotypes. It also highlights the importance of integration of molecular pathology signatures into clinical algorithms to optimize the future use of existing medications and inform the development of new drugs for refractory patients.

Unsupervised clustering defines treatment response diversity. Next, we used unsupervised analyses to explore the relationship of multiple genes/pathways with response to treatment. First, we applied principal component analysis (PCA) to identify underlying subgroup structures. PC1 and PC3 correlated with inflammatory cell infiltration in synovial biopsies, while they also associated with histological pathotypes primarily separating the lympho-myeloid and fibroid pathotypes (Extended Data Fig. 2a,b).
Unsupervised Monte Carlo consensus clustering (M3C) 11 showed 71% of rituximab responders (n = 24) in cluster 1 compared with only 29% (n = 10) in cluster 2 (P = 0.0004; Fig. 2a). Genes relevant for B cell biology were significantly higher in cluster 1 in patients treated with rituximab (Extended Data Figs. 2e and 3a). Cluster 1 was also linked with significant upregulation of the B cell gene module S136 from weighted gene correlation network analysis (WGCNA) 6 , together with upregulation of the proinflammatory M1 macrophage module S39 and downregulation of the fibroblast module S115 (Extended Data Fig. 2e).
Clustering of patients treated with tocilizumab was less distinctive, with 46% of responders (n = 21) in cluster 1 and 54% in cluster 2 (n = 25) (Fig. 2b). However, cluster 1 was significantly associated with IL-6 pathway genes (Extended Data Figs. 2f and 3b), together with upregulation of B cell and M1 macrophage modules and downregulation of fibroblast modules. In keeping with the increase in immune cell-related modules in cluster 1 for both treatments, semiquantitative IHC scores for synovial immune cells were significantly higher in cluster 1 (Extended Data Fig. 2g), indicating that immune cell infiltration is linked to gene expression in cluster 1, as inferred by the loss of significance when adjusting differentially expressed gene (DEG) analysis between consensus clusters 1 and 2 for immune cell content using PC1 as a covariate (Extended Data Fig. 3c,d). The strong correlation of PC1 with histology markers and immune cell-related genes (Extended Data Fig. 3e) is probably linked to this effect. Molecular signatures of treatment response. Next, we performed DEG analysis to identify genes associated with treatment response on all patients who at any point in the trial had received rituximab or tocilizumab (as described in Methods and Supplementary  Fig. 1). A total of 6,625 genes were significantly different (false discovery rate (FDR) < 0.05) in rituximab responders compared with nonresponders ( Fig. 2c and Supplementary Data 1), and 85 for tocilizumab ( Fig. 2d and Supplementary Data 1). Genes upregulated in the synovial tissue of rituximab responders included members of the immunoglobulin (Ig) superfamily and leukocyte-related genes. Fig. 1 | Synovial histological markers at baseline associate with response to rituximab and tocilizumab. a, Classification into synovial pathotypes according to semiquantitative scores for CD3 + T cells, CD20 + B cells, CD68 + macrophages and CD138 + plasma cells, with representative examples from patients classified as lymphomyeloid (CD20≥2 and/or CD138≥2), diffuse-myeloid (CD68SL≥2, and CD20/CD138<2) or fibroid/pauci-immune (CD68SL/CD20/CD138<2). Right, 16-week CDAI 50% response in patients stratified by pathotype (n = 152). Bar plots showing the proportion of CDAI 50% responders for rituximab (in blue) and tocilizumab (in yellow) within each pathotype, with corresponding exact numbers. Fisher's test, exact P values for P < 0.05. b, Approach to in silico deconvolution of synovial tissue using MCP-counter. c, MCP-counter scores for each cell type compared among CDAI 50% responders (R) and nonresponders (NR). Bar plots indicate nominal log 10 P values for tocilizumab and -log 10 P values for rituximab (two-sided Mann-Whitney test); dashed lines correspond to P = 0.05. Boxplots (right) show median and first and third quartiles, whiskers extending to the highest and lowest values. d-f, 16-week CDAI 50% response in patients stratified into B and T cell poor/rich (d) and macrophage/mDC poor/rich (e) according to median MCP-counter scores for individual cells (rich if above median, poor if below), or by combining B cell and macrophage/mDC scores from d,e (f). Exact P values shown when <0.05, two-sided Fisher's test comparing the proportions of responders to rituximab (in blue) and tocilizumab (in yellow). g-i, Longitudinal disease activity scores (CDAI), shown as mean ± s.d., for each month from baseline to 16 weeks for patients randomized to rituximab (in blue) or tocilizumab (in yellow) and classified as B and T cell poor/rich (g), macrophage/mDC poor/rich (h) and combined B cell/ macrophage poor/rich (i). Comparison of CDAI between the two medications at individual time points by two-sided Mann-Whitney test, exact P values for <0.05 (adjustment for multiple comparisons by FDR). P values for the drug × time interaction term (two-way repeated-measures analysis of covariance) are shown when <0.05. c-i, n = 133 patients with baseline RNA-seq. NK, natural killer cells. mDC, myeloid dendritic cells.  T cell rich T cell poor   Nonresponse to rituximab, on the other hand, was associated with complement genes, bone morphogenic proteins, fibroblast-related genes and several Hox genes. Interestingly, lymphocyte and Ig genes were also upregulated in the synovial tissue of tocilizumab responders. Both nonresponder groups showed upregulation of extracellular matrix genes, including integrin-binding sialoprotein, aggrecan and collagen, and genes linked to tissue remodeling, cell infiltration and cell-cell interaction. Following adjustment for immune cell infiltration by PC1, DEGs for rituximab remained significant and, in the case of tocilizumab, the number of identified DEGs increased (Extended Data Fig. 3f for rituximab and Extended Data Fig. 3g for tocilizumab; Supplementary Data1), suggesting that DEG analysis provides an additional dimension to the inflammatory cell infiltrate alone that differentiates responders from nonresponders. Of note, inclusion of covariates such as age, gender and ethnicity was not associated with major differences in the statistical significance of DEGs (Supplementary Data 1).
To investigate the functional role of the above genes, we applied quantitative set analysis for gene expression (QuSAGE) modular analysis 12 using blood-and synovium-specific WGCNA modules (Fig. 2e,f) 6,13 . Antigen presentation, T and B cell-related modules and interferon signaling were significantly increased in rituximab responders, while Hox gene and fibroblast modules were increased in rituximab nonresponders (Fig. 2e).
Myeloid cell cytokine, peroxisome proliferator-activated receptor (PPAR) and metabolic pathways were upregulated in tocilizumab responders (Fig. 2f). Although none of the modules was significantly modulated in nonresponders to tocilizumab, fibroblast modules were also detected in nonresponders to tocilizumab, suggesting the possible existence of a shared treatment-resistant signature.
Refractory disease is linked to a stromal/fibroblast signature. To further explore the hypothesis of a common refractory signature following treatment switch at 16 weeks ( Supplementary Fig. 1), we compared patients in whom both rituximab and tocilizumab failed to induce response (multidrug resistant/refractory, n = 40 for histology, n = 32 for RNA-seq) with (1) patients who responded exclusively to rituximab after tocilizumab failure (pro-rituximab, n = 11 for histology and n = 9 for RNA-seq) and (2) patients who responded exclusively to tocilizumab after rituximab failure (pro-tocilizumab, n = 13 for histology and n = 12 for RNA-seq) (Fig. 3a). We identified 1,980 genes upregulated in both pro-rituximab and pro-tocilizumab patients, 175 exclusive to the pro-rituximab group and 306 exclusive to the pro-tocilizumab ( Fig. 3b and Supplementary Data 2). Among genes upregulated in responders to both medications were lymphoid, myeloid and many cytokine genes (Fig. 3c,d). Chemokines and lymphocyte genes were upregulated in pro-rituximab patients, while lymphocyte and myeloid lineage genes were upregulated in pro-tocilizumab.
Modular analysis showed antigen presentation and dendritic, macrophage and plasma cell infiltration modules upregulated in responders to both biologics (Fig. 3e). Similarly, the CD8 and Tph T cell module was upregulated in each drug response group, with greater change for the rituximab responder group (proximity to pro-rituximab axis), while Toll-like receptor signaling and macrophage chemokine and cytokine signaling were significantly increased in pro-rituximab patients only. Modules for T cells, plasma cells and the TNF receptor superfamily gene were upregulated in pro-tocilizumab patients only (Fig. 3e).
Notably, 1,277 significant genes were unique to multidrugresistant/refractory patients (Fig. 3c,d and Supplementary Data 2), including fibroblast and extracellular matrix-encoding genes such as fibroblast growth factor (FGF), homeobox (HOX) and NOTCH family genes, together with multiple cell-adhesion-molecule-and collagen-encoding genes ( Fig. 3e and Supplementary Data 3).
In line with molecular signatures, baseline histological scores for CD3 + T cells and CD79a + B-cells and CD138 + plasma cells were significantly lower in refractory patients (Fig. 3f). Additionally, in silico deconvolution showed significantly lower levels of CD8 + T cells, monocytes and mDCs and a trend towards increase in endothelial cells, neutrophils and fibroblasts in refractory patients (Fig. 3g).
To further characterize the association of synovial fibroblast genes with multidrug resistance, we complemented MCP-counter deconvolution by examining enrichment in synovium-specific fibroblast gene modules derived from RA synovial single-cell RNA-seq 14 . As shown in Fig. 3h, the signature for HLA-DRA high sublining fibroblasts (SC-F2), a proinflammatory subset associated with leukocyte-rich synovial infiltration in RA, was significantly higher in responders (P = 0.027) as opposed to CD34 + sublining fibroblasts (SC-F1) and, in particular, to the newly described DKK3 + sublining fibroblasts (SC-F3), both increased in refractory patients (P = 0.036 and 0.00055, respectively).
For orthogonal confirmation of these findings at the protein level, we used multiplex immunofluorescence to detect DKK3 + fibroblasts in the synovial lining and sublining of refractory patients ( Fig. 3i and Extended Data Fig. 4).
Together, these results show that baseline histological and molecular signatures are associated with response to individual drugs, while nonresponse to multiple biologics is linked to a specific pretreatment signature associated with fibroblasts.
Digital spatial profiling of refractory RA. Because immune and stromal cells are known to exhibit positional identity relevant to the pathogenesis of RA 15 , we used digital spatial profiling (DSP)  Supplementary Fig. 1): patients responding to rituximab (RTX) following tocilizumab (TOC) failure (pro-rituximab, blue), patients responding to tocilizumab following rituximab failure (pro-tocilizumab, yellow) and patients in whom both drugs failed sequentially (refractory, red). Numbers in brackets denote patients with available RNA-seq. b, Venn diagram showing the overlap of DEGs between patients classified as in a. c,d, Three-way DEG analysis on baseline synovial biopsies of patients classified as in a, with side (c) and top view (d). Significant differences in pro-rituximab (blue), pro-tocilizumab (yellow) and refractory (red) patients and significant genes overlapping in pro-rituximab and pro-tocilizumab patients (green) are color coded. Significance was internally estimated by the volcano3D package combining significance (q < 0.05) from both LRT and pairwise Wald test via DESeq2. e, Three-way QuSAGE radial plot showing differential WGCNA module expression in patients classified as above. f, Histological semiquantitative scores for immune cells in refractory patients (n = 40) and responders to one of any two medications (n = 24). Boxplots showing median and first and third quartiles. Two-way Mann-Whitney test, exact P values FDR adjusted for multiple comparisons. g, Deconvolution of immune cells using MCP-counter in patients classified as refractory or responders as in a. Boxplots showing median and first and third quartiles, dot-plots showing individual patients. Two-way Mann-Whitney test, exact P values FDR adjusted for multiple comparisons. h, Fibroblast single-cell subset enrichment scores in refractory patients (n = 32) or responders to either rituximab or tocilizumab (n = 21), as in a. Boxplots showing median and first and third quartiles, whiskers extending to the highest and lowest values. Exact P values are shown, two-sided Mann-Whitney test. i, Multiplex immunofluorescence in refractory and responder patients; nuclear staining (blue), CD45 (red), CD90 (green), DKK3 (yellow) (all top) and DKK3 single staining (yellow, bottom). *, DKK3 + CD45 + lymphocytes; arrowheads, DKK3 + CD90 + fibroblasts. A larger overview and individual stainings are provided in Extended Data Fig. 4. Representative images out of a total of three refractory and three responders. Scale bars, 50 μm. NS, not significant.
to characterize the spatial positioning of cell signatures in association with treatment response/resistance. We employed GeoMx DSP (NanoString), which uses a set of protein lineage markers to define regions of interest (ROIs) that undergo whole-transcriptomic spatial RNA expression (Fig. 4a). First, we compared gene expression in responders and refractory patients across all ROIs: lining/ superficial sublining, deep sublining and lymphoid aggregates (Fig. 4b). Consistent with the above bulk RNA-seq modules and protein expression, multiple genes related to the DKK3 + fibroblast subset (PRELP, OGN, CAM1KD) were significantly higher in refractory patients (Fig. 4c). When looking at individual ROIs, we found specific genes differentially expressed in responders versus refractory patients in each synovial region. For example, the gene encoding for the fibroblast marker FAP 16 was significantly upregulated in the deep sublining of refractory patients, the gene for the osteoclast marker RANK (TNFRSR11A) was significantly higher in 1.2 S102 Myofibroblast S114 Neutrophil degranulation S116 Fibroblast_1_CD55 + S185 Hox genes S191 Chondrocyte differentiation S194 Platelets S200 Fibroblast_2a_THY1 + S119 Mast_cells, neutrophil degranulation S120 Plasma_cells S126 Macrophage S13 Plasma_cells, XBP1 signaling S139 CD8 and Tph T_cells S156 Dendritic_cells, MHC class I antigen presentation S166 Macrophage, Toll-like receptor (TLR) signaling

Pre-and post-treatment histopathological and molecular analyses.
To explore the longitudinal effects of each drug on synovial immune cell infiltration and gene expression, we compared paired synovial samples at baseline and 16 weeks (rituximab, n = 41 for histology and n = 29 for RNA-seq; tocilizumab, n = 24 and n = 15, respectively) (Extended Data Fig. 5a and Extended Data Table 1). First, by histology, we showed a significant reduction in synovial CD20 + total B cells, CD79a + B cells and CD138 + plasma cells in patients treated with rituximab, in line with the rituximab mechanism of action targeting CD20 + B cells (Fig. 5a). Conversely, patients treated with tocilizumab showed a significant reduction in CD68 + sublining macrophages but not B-cells (Fig. 5a). Analysis of covariance showed a significantly higher reduction of CD20 + and CD79a + B cells in patients treated with rituximab, and a significantly higher reduction in CD68 + sublining macrophages in those treated with tocilizumab (Extended Data Fig. 5b).
Similar results were obtained when comparing MCP-counter immune cell signatures. Namely, patients treated with rituximab showed a significant reduction in B cells, T cells and monocytes/ macrophages while those treated with tocilizumab showed a significant reduction in monocytes/macrophages and T cells, but also in neutrophils and mDCs and, interestingly, an increase in fibroblast signature (Extended Data Fig. 5c). This suggests that both biologics have an effect on immune cells but that tocilizumab can potentially also affect stromal cells.
To further dissect these longitudinal molecular signatures, we developed an R package to fit negative binomial mixed-effects models at the individual gene level (general linear mixed-effects model (glmmSeq)), because mainstream RNA-seq analysis tools are unable to fit mixed-effects linear models (Methods).
Using glmmSeq to compare gene expression over time in paired synovial biopsies, 7,316 genes were significantly up-or downregulated by both drugs while 345 were differentially affected by either drug based on significance (FDR < 0.05) of the interaction term time × medication ( Fig. 5b and Supplementary Data 4a). Of note, MS4A1 (encoding CD20), PAX5 and BLK were significantly downregulated in response to rituximab, consistent with B cell depletion mechanism and histology results (Fig. 5c), while tocilizumab induced a reduction in IL-6-related transcripts, also consistent with the IL6 acting mechanism of tocilizumab, not the CD20 mechanism of rituximab.
When patients were stratified according to response, a significant reduction in CD138 and CD79a plasmablasts/plasma cells was observed only in rituximab responders, while a significant reduction in CD68SL macrophages was observed only in responders to  tocilizumab (Extended Data Fig. 5d), indicating that reduction in B and synovial plasma cells and macrophages is associated with response to rituximab and tocilizumab, respectively. The mixed-effects model allowed us to further examine various changes in gene expression following therapy between responders and nonresponders to each drug (Extended Data Fig. 5e-h). Rituximab had a general effect on 1,796 genes, with 349 showing significant (FDR < 0.05) differential expression change over time between responders and nonresponders (Extended Data Fig. 5e and Supplementary Data 4b). Rituximab responders showed a greater decrease in SAA1 and SAA2 (serum amyloid proteins 1 and 2), as well as greater decreases in Ig chain genes IGHV3-64D and IGKV1-13, suggesting that a drop in antibody-secreting B cells is associated with response to rituximab (Extended Data Fig. 5g). Chemokine-encoding CXCL11, the citrullination enzyme encoding PADI2 (peptidyl-arginine-deiminase2) and the key Th17 and mucosal-associated invariant T (MAIT) cell transcriptional regulator RORgamma (RORC) gene were also modulated in rituximab responders (Extended Data Fig. 5g).
Tocilizumab treatment resulted in modulation of 1,609 genes, with an additional 136 showing differential change in gene expression between responders and nonresponders (Extended Data Fig. 5f and Supplementary Data 4c). Reduction in pro-lymphoid follicle development genes encoding for lymphotoxin-A (LTA), complement receptor-2 (CR2), lymphoid-tissue-resident dendritic cell marker XCR1 17 and prolactin (CLEC17A), expressed on proliferating germinal center B cells 18 , augured response to tocilizumab (Extended Data Fig. 5h).
To further investigate pathway modulation induced by treatment, genes identified in the longitudinal mixed-effects model analysis were analyzed for Gene Ontology (GO)/pathway enrichment. Rituximab treatment induced significant downregulation of B cell receptor pathways, as well as bone resorption and remodeling pathways (Fig. 5d). When stratifying patients according to response, responders to rituximab showed significant downregulation of T cell receptor complex, lymphocyte chemotaxis and migration, chemokines and IL-1-related pathways (Fig. 5e), suggesting that rituximab response is linked to additional immunomodulation in addition to reduction in B cell-related pathways. Responders to tocilizumab showed significant decrease in humoral immune response, Ig, B cell and complement activation (Fig. 5f), in line with the known effect of IL-6 on B cell growth/differentiation.
In summary, longitudinal analyses of matched pre-and post-treatment biopsies indicate that specific biological changes are associated with response to individual treatments.
Machine learning models predict drug response and multidrug resistance. To establish the ability of synovial tissue gene expression in prediction of treatment response/resistance, we developed machine learning predictive models with the dataset partitioned for training and testing using ten-by-tenfold nested cross-validation, as detailed in Methods and schematically in Fig. 6a. Supplementary Table 1 shows the performance of models used to predict (1) rituximab response, (2) tocilizumab response and (3) refractory state. Final models (Supplementary Table 2) were trained on the entire dataset to extract variable importance (Fig. 6a, bottom right and Extended Data Table 2).
As shown in Fig. 6b, the optimal predictive models included gene elastic net regression for rituximab and tocilizumab response, with 40 genes (AUC = 0.744) and 39 genes (AUC = 0.681), respectively, and gradient-boosted machine (GBM) for refractory state, with 53 genes (AUC = 0.686) (Extended Data Fig. 6a). Notably, no clinical features were selected by the final models (Extended Data Table 2) and, in comparison with RNA-seq, predictive models built using clinical and histology parameters alone performed quite poorly (Extended Data Fig. 6b). AUC values in the omitted inner cross-validation folds were consistent with AUC results in the true test folds. Multiple genes were shared across models, with only 85 required to build all three prediction models and 32 shared between at least one model (Extended Data Fig. 6a and Extended Data Table 2). Each model selected multiple genes of biological relevance to synovial tissue inflammatory and repair responses, as well as to bone and cartilage biology. Key prediction genes shared between all three models included genes encoding for: XCR1, a marker of DC1 migratory DCs; 17 chemokine CXCL14; acute phase reactant SAA2; and IGHV7-4-1, probably reflecting tissue-resident plasma cells. The refractory state model, which contained the largest number of unique genes, included several linked to the fibroid pathotype such as TNFRSF11B, which encodes the osteoclast negative regulator osteoprotegerin and the chondrocyte adhesion mediator CHAD (chondroadherin), but also the citrullination enzyme PAD4 encoding gene PADI4 19 , consistent with a role of persistent tissue destruction and remodeling in the refractory RA state.
Although the original clinical trial was not powered for these types of analyses, and larger cohorts will be required to further validate the models and improve their predictive power, these results show that predictive models can harness molecular information from synovial biopsies at baseline before treatment and thus are of potential clinical utility for prediction of response to therapy.

discussion
Our study provides an in-depth molecular and histological profiling of joint tissue from a biopsy-driven, randomized clinical trial in RA (R4RA) 9 , affording insights into the cellular and molecular pathways underpinning the diverse treatment response to two commonly used targeted biologic therapies directed against B cells (rituximab) and the IL-6 receptor (tocilizumab). Using both conventional histology and in silico deconvolution, we observed that lymphoid cells were associated with response to rituximab while myeloid cells were associated with response to tocilizumab. Although the importance of synovial macrophages as a predictor of response to anti-TNF has been described 20,21 , here we report that within the synovial B cell-poor group, an enhanced response to tocilizumab is associated with the presence of myeloid cells. In regard to rituximab, previous small observational studies reported pretreatment synovial CD79a + B cells 22 and synovial molecular signatures 23 as potential response predictors.
In this larger, randomized clinical trial cohort we confirmed some of these findings and, combining histological assessment with advanced molecular analyses, we identified genes and pathways linked to the cognate drug targets in association with response. For rituximab these were B cell genes, Igs, chemokines and leukocyte genes. Response to tocilizumab was associated with IL-6 pathway genes, but also with lymphocyte and immunoglobulin genes, which is not surprising since IL-6 is a well-known B cell growth factor 24 .
Modular analyses demonstrated how genes increased in rituximab responders were functionally related to antigen presentation and lymphocyte activation, together with interferon signaling genes, in keeping with previous reports linking increased type I interferon with response to rituximab 25,26 . In tocilizumab responders, in line with the prominent role of myeloid cells identified by histopathology, the myeloid cell cytokine module was upregulated together with PPAR signaling and metabolic pathways.
In contrast, nonresponse to both drugs was defined by >1,000 genes and several shared pathways, including Hox, FGFs and ECM genes/modules. Critically, this shared nonresponse signature is linked to the fibroid pauci-immune pathotype characterized by scanty immune cell infiltrate with prevalence of stromal cells 6,7 , which we also reported as being associated with poor response to synthetic disease-modifying antirheumatic drugs (DMARDs) 27 and TNF inhibitors 28 . This supports the concept that the fibroid pauci-immune phenotype represents a refractory endotype, since   Fig. 6 | Predictive models using nested ten-by-ten-fold cross-validation for response to rituximab and tocilizumab. a, Machine learning pipeline utilized to predict CDAI 50% response to rituximab and/or tocilizumab using gene expression, clinical data and histological data as features (n = 133). Data processing (1) involved selection of protein-coding genes with the highest variance and removal of highly correlated genes. Data were split into ten inner and ten outer folds for building machine learning models (2). In models built using gene expression, RFE or univariate filtering was used to select the most important/predictive features for each model. Each model was evaluated on both the test set and the set omitted during cross-validation (3). Average tuned parameters from the outer folds were used to fit to the whole dataset to determine the importance of features selected for each model (4). b, Grid of plots showing optimal predictive models for different treatments (left, glmnet rituximab response prediction; middle, glmnet tocilizumab response prediction; right, GBM refractory response prediction) using gene expression and baseline clinical parameters as features. From top to bottom, plots show ROC curves for the best model on the test dataset (from outer fold), ROC curves on the omitted dataset (from inner fold) and variable importance when fit to the whole dataset.
in R4RA double nonresponder patients conventional synthetic DMARDs (csDMARDs) and TNF inhibitors had already failed, as per trial entry criteria, thus displaying resistance to three biologic therapies targeting distinct immunological pathways (TNF, CD20 + B cells and IL6R) and meeting the definition of multidrug resistant/ refractory RA 2 . Because recent studies using single-cell RNA-seq (scRNA-seq) on RA synovium identified specific fibroblast subsets with critical roles in RA pathogenesis 13,14,16,29 , we applied modular analysis based on scRNA-seq fibroblast subsets 14 and orthogonal validation by multiplex immunofluorescence, and identified an association between DKK3 + fibroblasts and refractoriness. DKK3 encodes Dickkopf3, a negative regulator of beta-catenin that has been shown to promote aggressive behaviour in cancer-associated fibroblasts 30 , although the exact function of DKK3 + fibroblasts in RA remains to be established 31 . Furthermore, DSP revealed specific upregulation of genes in the sublining of refractory patients encoding for the fibroblast marker FAP, which has been linked with RA pathogenesis 16 , while other markers consistently modulated across all regions included CCL13 encoding for the monocyte-attracting chemokine MCP-4, which has been shown to activate synovial fibroblasts 32 . Hence, stromal cells-and, in particular DKK3 + fibroblast genes-may be a new drug target that helps overcome the complex problem of refractoriness in RA.
Genes linked with response included increased expression in the synovial aggregates encoding for CD24, a B cell marker associated with response to biologic treatment 33 . Also consistent with the capacity of biologic therapies to halt structural damage progression 34 , we identified increased expression of RANK transcripts (TNFRSF11A) in the lining/superficial sublining of responders, in line with the reported presence of osteoclast precursors in inflamed synovia 35 . Overall, these results suggest that the spatial organization of immune infiltrates is highly relevant for determinion of treatment response/resistance, although additional work is needed to dissect the contribution of specific markers in individual synovial regions and their association with therapy response/resistance.
The longitudinal analysis of matched pre-/post-treatment synovial biopsies enabled us to investigate drug effects on synovial pathology and gene expression. Rituximab reduced synovial CD20 + B cells in both responders and nonresponders but, notably, a clinically relevant response required broader and deeper impacts on differentiated CD79a + plasmablasts and CD138 + plasma cells over and above CD20 + B cell depletion. These observations are in agreement with previous observational studies showing changes in plasma cells between responders and nonresponders to rituximab 36 . However, another study showed variable depletion of B cells and plasma cells and an unclear association with treatment response 37 , although the small sample size and biopsy analysis at different time points (4, 8, 12 and 16 weeks) make it difficult to draw conclusions. In our study, because repeated synovial biopsy was performed at 16 weeks following one cycle of rituximab (2× 1-g infusion), it is plausible that a repeated biopsy at the time of the second infusion at 6 months (not available in sufficient numbers to be informative) could have detected a wider/deeper effect on B cell lineages linked to clinical response.
Following tocilizumab therapy, nonresponders were characterized by a failure to reduce sublining macrophages, which is consistent with previous literature indicating macrophages as markers of treatment response 20 .
To assess the relationship of gene expression changes and treatment response, we developed a pipeline for mixed-model analysis of RNA-seq data in repeat biopsies. This revealed patterns of change in gene expression not detectable by standard analytical pipelines, while interaction analysis allowed us to identify genes that were affected by each drug specifically. Biological differences in synovial gene expression following treatment with rituximab or tocilizumab were consistent with the cognate treatment targets: B cell depletion and IL-6 receptor blockade, but also revealed unexpected differences such as differential changes in metalloproteinases. Interaction analysis showed that rituximab responders had a greater decrease in genes encoding for serum amyloid proteins, Ig chains, the citrullination enzyme PADI2 and transcriptional regulator RORgamma. In tocilizumab responders, a greater reduction in pro-lymphoid follicle development genes was observed in keeping with the important role of IL-6 as a B cell growth factor driving in situ ectopic lymphoid structure development within inflamed tissues 38 .
In translational terms, the importance of molecular studies is measured by their ability to enhance disease understanding, but also on their clinical impact 39 . Thus, to determine the predictive value of deep molecular characterization in foretelling treatment response, we applied a number of machine learning methodologies resulting in the selection of models effective at predicting treatment response as tested by nested cross-validation (AUC for rituximab = 0.744, tocilizumab = 0.681, refractoriness = 0.686). The purpose of testing multiple models was to determine whether nonlinear algorithms (least-squares support vector machine (SVM), penalized discriminant analysis (PDA), mixture discriminant analysis (MDA), tree-based and so on) could outperform elastic net linear regression. In practice, for prediction of response to each drug individually, elastic net regression performed best. However, a GBM model was superior in predicting refractoriness, consistent with the notion that biological heterogeneity underlying refractoriness to multiple drugs might require a nonlinear algorithm for optimal prediction. Of relevance to future clinical practice, gene expression models were clearly superior to those built using clinical and histological data alone. Seropositivity, which has been weakly associated with response to rituximab 40 was not included in the final models.
A known limitation of this study is the relatively modest sample size of the training data for each predictive model. However, the study was powered purely for the primary outcome of the original clinical trial and numbers reflect the difficulty in conducting biopsy-driven randomized trials in RA. With the current sample size and case/biomarker ratio, statistical theory suggests that we are likely to be only part of the way up the learning curve 41 . The AUCs derived in this study would still be relatively low for direct application to clinical use without further validation and improvement. Refining and reducing the number of genes in the models with validation in independent cohorts could improve the prediction models, as could incorporation of information from single-cell studies and enhanced deconvolution methods. Another limitation with regard to interpretation of the longitudinal results is that the second biopsy at 16 weeks was an optional procedure and, thus, at risk of selection bias. Although there were no major baseline differences, patients who underwent repeated biopsy had lower response rates, which is expected because responders would have been less likely to consent to a second biopsy.
In conclusion, this study provides insights from analysis of diseased tissue regarding the mechanisms driving treatment response heterogeneity in RA, and underscores the importance of integrating predictive molecular pathology signatures into clinical algorithms to optimize the usage of existing drugs. The identification of genes and cell types associated with multidrug resistance could aid the development of new drugs for refractory patients in whom current medications targeting classical immune pathways are not effective. We envisage that routine use of synovial biopsies could facilitate a patient-centered approach 5 to the management of RA, thus moving away from the current trial-and-error drug prescribing towards an emergent era in which selection of the optimal drug is based on synovial biopsy gene signatures.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of

Patients and intervention.
A total of 164 patients aged 18 years or over, fulfilling the 2010 American College of Rheumatology/European Alliance of Associations for Rheumatology (EULAR) classification criteria for RA and who were eligible for treatment with rituximab therapy according to UK NICE guidelines-that is, failing or intolerant to csDMARD therapy and at least one biologic therapy (excluding trial IMPs), were recruited when fulfilling the trial inclusion/exclusion criteria (for the full study protocol and baseline patient characteristics see Humby et al. 9 ). Briefly, patients underwent synovial biopsy of a clinically active joint at entry to the trial, performed according to the expertise of a local center as either an ultrasound-guided or arthroscopic procedure 42 Following synovial biopsy, patients were randomized to receive rituximab as two 1,000-mg intravenous infusions 2 weeks apart or intravenous tocilizumab at a dose of 8 mg kg −1 administered at 4-weekly intervals. Patients were followed up every 4 weeks throughout the 48-week trial treatment period, during which RA disease activity measurements and safety data were collected. An optional repeated synovial biopsy of the same joint sampled at baseline was performed at 16 weeks (Supplementary Tables 3  and 4). The study protocol has been published online (http://www.r4ra-nihr. whri.qmul.ac.uk/docs/r4ra_protocol_version_9_30.10.2017_clean.pdf) and was registered on the ISRCTN database (no. ISRCTN97443826) and with EudraCT (no. 2012-002535-28). Patient demographics are reported in the original trial publication 9 . All patients provided written informed consent. Participants did not receive any compensation, except for reimbursement of travel expenses. Response criteria and treatment switch. The primary endpoint was defined as CDAI ≥ 50% improvement from baseline at 16 weeks 9 . CDAI is calculated by totaling the number of tender joints (0-28), the number of swollen joints (0-28), patient global health assessment on a 0-10 visuoanalogic scale and the care provider global health assessment on a 0-10 visual analog scale.
As shown in Supplementary Fig. 1, CDAI 50% nonresponders at 16 weeks were switched to the alternative biologic agent and their response was assessed at 16 weeks following the switch, as determined by CDAI 50% improvement. Including crossover patients, a total of 108 patients were treated with rituximab and 117 with tocilizumab. Of those treated with rituximab, 43 were defined responders (40%) while 53 responded to tocilizumab (45%). Among all responders, 11 responded to rituximab following tocilizumab failure and were classified as exclusive responders to rituximab (pro-rituximab), while 13 responded to tocilizumab following rituximab failure and were thus classified as pro-tocilizumab. Patients in whom both drugs failed throughout the study were classified as multidrug resistant/refractory (n = 40).

RNA-seq and molecular classification/analysis.
A minimum of six synovial samples per patient were immediately immersed in RNA-Later and RNA was extracted from tissue using one of two protocols: phenol/chloroform isolation and Zymo Direct-zol RNA MicroPrep-Total RNA/miRNA Extraction kit. In both methods, tissue was lysed in Trizol solution using a LabGen125 homogenizer. Briefly, for the phenol/chloroform extraction method, 1-10 mg of tissue was lysed and then sheared using a 21 G needle. The tissue lysate was then mixed vigorously with chloroform before centrifugation. The aqueous phase was removed and mixed with ice-cold isopropanol for 30 min. Following further centrifugation, the RNA pellet was washed in 70% ethanol before air-drying and resuspension in RNAse-free water. Samples extracted using Zymo Direct-zol Miniprep kits were processed as per the manufacturer's instructions. Briefly, 1-10 mg of tissue lysate was run through the Zymo-Spin IC column. Columns were then washed using the appropriate kit wash buffers before RNA was eluted and resuspended in RNAse-free water. Quality control was carried out by quantifying samples via spectrophotometer readings on a Nanodrop ND2000C. RNA integrity was measured using Pico-chip technology on an Agilent 2100 Bioanalyzer to determine RNA integrity number. A total of 214 synovial tissue samples were available for RNA extraction and were subsequently sent for RNA-seq to Genewiz. RNA-seq libraries were prepared using the NEBNext Ultra RNA Library Prep kit for Illumina, following the manufacturer's (NEB) instructions. Briefly, messenger RNAs were initially enriched with Oligo d(T) beads followed by limited PCR cycles. The sequencing library was validated on an Agilent TapeStation (Agilent Technologies) and quantified using a Qubit 2.0 Fluorometer (Invitrogen), as well as by quantitative PCR (KAPA Biosystems). The sequencing libraries were clustered on Illumina flowcells. Sequencing was performed on an Illumina HiSeq instrument according to the manufacturer's instruction. Samples were sequenced using a 2 × 150-base-pair (BP) paired-end configuration.
RNA-seq data processing. A total of 214 paired-end RNA-seq samples from 50 million reads of 150-bp length were trimmed to remove the Illumina adapters using bbduk from the BBMap package v.37.93, with default parameters. Transcripts were then quantified using Salmon 46 v.0.13.1 and an index generated from the Gencode release 29 transcriptome following the standard operating procedure. Tximport v.1.13.10 was used to aggregate transcript-level expression data to genes, then counts were subjected to variance-stabilizing transformation (VST) using the DESeq2 v.1.25.9 package 47 . Following RNA-seq quality control, 36 samples were excluded due to poor mapping or RNA quality. Using unsupervised PCA and plotting the first five eigenvectors in pairs, one outlier was identified and removed from further analysis. Thus RNA-seq data from 133 patients were available for subsequent analysis at baseline, and from 44 patients for the follow-up time point. Baseline characteristics of patients with available RNA-seq are shown in Supplementary Table 5. The first six PCs did not associate with demographics, treatment and its associated response or clinical disease features such as disease activity or anticyclic citrullinated protein antibody status (Extended Data Fig. 2a,c,d).
Starting with length-scaled transcripts per million (TPM) counts derived using the R package tximport, Limma voom was used for normalization of data and calculation of weights for linear modeling 48 .
Cluster analysis. For cluster analysis, after removal of low-expressed genes, VST data were filtered using a coefficient of variation cutoff of >0.075 to select the 22,256 (of 56,809) most variable genes . These genes were used for cluster analysis of all baseline patients (n = 133) using the M3C algorithm 11 with partitioning around medoids clustering and 1,000 iterations. The lowest penalized cluster stability Index was used to select the number of clusters. After cluster assignment, patients were split into treatment groups using Pearson's distance metric and complete linkage method, and plotted using the ComplexHeatmap package (v.2.2.0) in R. An χ 2 test was applied to test significance between clusters and response to treatment based on the trial primary outcome measure (CDAI 50%) and, additionally, EULAR C-reactive protein (CRP) response (EULAR response) as another commonly used criterion.
Differential expression and modular analysis of RNA-seq data. Samples from all patients treated with either rituximab or tocilizumab throughout the trial were included in DEG analysis. This also comprised nonresponders who, as per trial protocol, were switched to the alternative medication at week 16, as shown in Supplementary Fig. 1. Neither responders nor nonresponders showed any significant differences in their baseline characteristics, including histological and molecular B cell status, gender or disease duration (Supplementary Table 6). Low-expressed genes were excluded from analysis, with the remaining 30,841 used for DEG analysis. This was based on negative binomial distribution via regression models of normalized count data using DESeq2, and a Wald test to compare variation between treatment response groups in synovium RNA-seq samples. Wald test-derived P values were FDR adjusted using Storey's q-value, with a cutoff of q < 0.05 used to define significantly DEGs (Supplementary Data 1). Distributions of DEGs are illustrated in volcano plots, and DESeq2 outputs were used for further modular analysis with the Bioconductor package QuSAGE v.2.10.0. Gene modules from Li et al. 12 and WGCNA modules were selected for gene set enrichment.
Deconvolution. MCP-counter 10 was used to deconvolute synovial RNA-seq, with the package Immunedeconv. Following deconvolution, patients were classified into rich/poor according to the median value of the individual cell type (for example, B cell rich if above the median value of MCP B cells, poor if below). For the enrichment of four fibroblast subtypes (SC-F1: CD34 + sublining, SC-F2: HLA + sublining, SC-F3: DKK3 + sublining and SC-F4: CD55 + lining), we used average expression of gene signatures obtained from differential gene expression analysis and known markers previously described by scRNA-seq 14 . Module scores for each subtype were calculated using the AddModuleScore function in the R package Seurat. The top five differentially expressed genes were considered subtype-specific gene sets and did not have genes in common. Wilcoxon testing was used for statistical assessment of module scores when comparing responders and nonresponders.
Crossover analysis of patients who underwent treatment switch. The drug-crossover analysis was performed on baseline RNA-seq samples of patients who underwent treatment switch (Fig. 3a). RNA-seq counts of protein-coding genes (n = 19,508) were used to perform a likelihood ratio test (LRT) that was calculated in comparison to a reduced model with the DESeq2 R package. Three-dimensional volcano plots and radial plots were generated using the volcano3D (v.1.0.3) package in R (Fig. 3c-e). QuSAGE was applied using WGCNA-derived gene modules, and radial plots were created using the volcano3D package with a P value significance threshold of <0.05 (Fig. 3e).
Multiplex immunofluorescence. Immunofluorescence staining was performed on 3-µm, formalin-fixed, paraffin-embedded human sections obtained from synovial tissues of patients with RA. Tissue sections were deparaffinized in sequential changes of xylene and ethanol chambers before washing and placing in a preheated target retrieval solution (pH 6.0; Dako, no. S1699) in a pressure cooker for 15 min. Tissue sections were allowed to cool at room temperature (RT) before washing in Tris-buffered saline (TBS). Endogenous peroxidase and biotin activity were blocked with peroxidase (Dako, no. S2023) for 10 min at RT.
DAPI (Thermofisher) nuclear counterstaining was applied for 10 min at RT and slides were then mounted with ProLong Gold Antifade reagent (Thermofisher).
For the NanoString GeoMx DSP WTA assay, slides were prepared following the automated Leica Bond RNA Slide Preparation Protocol (NanoString, no. MAN-10131-03). In situ hybridizations with the GeoMx Whole Transcriptome Atlas Panel (WTA, 18,677 genes) at 4-nM final concentration were done in Buffer R (NanoString). Morphology markers were prepared for four slides concurrently using Syto13 (DNA), CD20, CD3 and CD68 in Buffer W for a total volume of 125 μl per slide. Slides incubated with 125 μl of morphology marker solution at RT for 1 h, then washed in SSC and loaded onto the NanoString DSP instrument.
On the DSP instrument each slide was scanned with a ×20 objective at scan parameters 60 ms FITC/525 nm, 200 ms Cy3/568 nm, 250 ms Texas Red/615 nm and 300 ms Cy5/666 nm.
The resulting immunofluorescent images were used to select six freeform polygon-shaped ROIs containing approximately 200 nuclei in CD68 + synovial tissue lining and superficial sublining, CD20 -CD3sublining and CD20 + CD3 + lymphocyte aggregates.
After approval of ROIs, GeoMx DSP photocleaved the ultraviolet (UV)-cleavable barcoded linker of bound RNA probes and collected individual segmented areas into separate wells in a 96-well collection plate.
The dataset included 72 ROIs from 12 patients (four refractory and eight responder) across the three ROI types. An NTC water well was used for quality control checks. DSP analysis. GeoMx WTA sequencing reads from NovaSeq6000 were compiled into FASTQ files corresponding to each ROI. FASTQ files were then converted to digital count conversion files using the NanoString GeoMx NGS DnD Pipeline. Out of 18,677 genes, 17,065 exceeded the lower level of quantitation (LOQ) in >10% of ROIs; genes that did not exceed LOQ were excluded from the analysis. For normalization, counts were divided by sample-specific size factors determined by the median ratio of gene counts relative to geometric mean per gene. The DESeq2 R package was used for this preprocessing step. Differential expression analysis. We conducted differential expression analysis to compare responders and refractory patients using DESeq2 50 . This analysis was done for all ROIs simultaneously (responders, n = 48; refractory, n = 24) but separately for each location in the synovial layer: CD68 + lining/superficial sublining (responders, n = 17; refractory, n = 8); CD20 -CD3deep sublining (responders, n = 21, refractory n = 12); and CD3 + CD20 + lymphoid aggregates (responders, n = 10; refractory, n = 4). Since samples were collected from different locations, in the analysis of all samples we included location as a covariate (~location + response) to eliminate its influence on gene expression. The qvalue R package implementing Storey's q-value method was used to correct for multiple testing effects, and a cutoff of q < 0.05 was used to define significantly DEGs.
Longitudinal mixed-effects model analysis. Longitudinal analysis of RNA-seq on paired synovial biopsies was performed by fitting a negative binomial distribution GLMM for each gene. Because the most widely used mainstream differential gene expression analysis tools-edgeR, DESeq2 47 and Limma voom 48 -are all unable to fit mixed-effects linear models, we developed the R package glmmSeq to fit negative binomial mixed-effects models at the individual gene level. glmmSeq uses the glmer function from the R package lme4 (v.1.1-25), with negative binomial family function from the MASS package (v.7.3-53). Models were fit using maximum-likelihood estimation by Laplace approximation and bound optimization by quadratic approximation. For analysis of the differential effects of the two trial medications over time, the following model was fitted for each gene individually: where Y ijg is the longitudinal raw count of gene g in individual I at timepoint j, α g is the dispersion parameter for each gene, o ij is an offset term scaled to the logarithm of the total library size for each sample, b gi are random effects between individual patients, and N and NB are the functions for normal distribution and negative binomial distribution, respectively. TPM counts were used as input, and only individuals with paired samples were included (88 samples, 44 individuals The R package glmmSeq is downloadable via CRAN and the source code is also available, from https://github.com/KatrionaGoldmann/glmmSeq. When compared against a Gaussian linear mixed-effects model on log count data, glmmSeq showed similar results with strong correlation between P values generated using either distribution ( Supplementary Fig. 2a-c). Q-Q plots suggested that the negative binomial mixed model showed greater power in identification of significant effects ( Supplementary Fig. 2d-f).
Building classifier models for prediction of rituximab and tocilizumab response and refractory status. Baseline gene expression and clinical and histological data were used as features for machine learning models built to predict CDAI 50% response to either rituximab or tocilizumab treatment at the primary endpoint (16 weeks) or refractory status, defined as nonresponse to both drugs at the secondary endpoint (post-treatment crossover, 24 weeks). An overview of the pipeline is shown in Fig. 6a.
Although the R4RA study was not powered for machine learning, information theory demonstrates that sparse models developed from large biomarker panels in which only a small percentage of biomarkers have nonzero effects can still demonstrate evidence of prediction with relatively modest sample size, although the small sample sizes in our study mean that the predictive models are likely to be only part of the way up the learning curve 41 . In the present study ~2% of the 1,500 biomarkers inputted into the modeling system have a nonzero effect. If the C-statistic of the optimal classifier is 0.83, a sample size equivalent to 0.05 events per variable (n = 83) would be required to learn a classifier that has expected information for discrimination equal to 25% of that obtained by the optimal classifier (equivalent to a C-statistic of 0.68) (see the online calculator, https://pmckeigue.shinyapps.io/sampsizeapp/) 41 . The model feature space was created using either clinical and histological parameters or clinical data with gene expression. Gene expression data underwent VST and were subset to protein-coding genes (using gencode gene annotation v.29) with the highest expression variance (top 10%). Highly correlated genes (r > 0.9) were removed using the findCorrelation function from the R package caret (v.6.0-86), leaving 1,438 genes. Clinical features included: baseline tender joint count (TJC), swollen joint count (SJC), age, gender, CDAI, erythrocyte sedimentation rate (ESR), CRP and disease activity score based on ESR and CRP (DAS28ESR and DAS28CRP, respectively). Histology features included CD3, CD68L, CD68SL, CD20 and CD138.
Following processing, data were split into 10 × 10 nested folds (Fig. 6a  (2)). For models using gene expression features, filtering was performed using either RFE or univariate filtering from the caret package v.6.0. The number of features selected was chosen to maximize accuracy from 25, 30, 50 or 100. Model hyperparameters were tuned by inner tenfold cross-validation, with model accuracy determined in separate outer cross-validation folds to give an unbiased estimate of model accuracy.
Seven machine learning methodologies from the caret package were used to create the classifier models: elastic net (glmnet), random forest (RF), least-squares support vector machine (SVM) with radial basis function kernel (svmRadial), least-squares SVM with polynomial kernel (svmPoly), GBM, MDA and PDA. Models that failed to converge during training were excluded from evaluation. The purpose of testing multiple models was to determine whether nonlinear decision boundaries-as used by SVM, MDA, PDA and tree-based prediction algorithms such as GBM-could outperform penalized linear regression.
To evaluate model performance, receiver operating characteristic (ROC) curves were built using the plotROC R package v.2.2.1 to determine prediction accuracy in the outer fold test data and samples omitted for the inner fold. AUC was calculated to determine prediction performance. Tuning parameters for the final model were finalized as the mean over all ten outer folds. The final best model for each classification was fit to the entire dataset, exported and feature importance ranked.

Statistical analysis.
For cross-sectional comparisons of continuous variables between two groups the Mann-Whitney U-test was used, whereas the Wilcoxon signed-rank test was used to assess the difference between groups with longitudinal paired data. More specific analyses of RNA-seq count data are detailed above in each relevant section. R v.4.0.0, or later, was used for all formal testing analyses.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

data availability
The datasets generated during and/or analyzed during the current study are available on an interactive web interface that allows direct data exploration (https:// r4ra.hpc.qmul.ac.uk/). A searchable interface is available to examine relationships between individual synovial gene transcript levels and histological and clinical parameters, and clinical response at 16 weeks. In addition, interactive versions of Figs. 3c and 5b and Extended Data Fig. 5e,f allow users to click on individual genes to see their expression and search for genes of interest. The website was constructed using R shiny server 1.5.16, with interactive plots generated with R plotly 4.9.3. The datasets can be downloaded from https://www.ebi.ac.uk/arrayexpress/ experiments/E-MTAB-11611. Other public datasets used for pathway analysis were sourced from the GO annotation (GOA) database (BiologicalProcess-EBI-UniProt-GOA (11 February 2020), CellularComponent-EBI-UniProt-GOA (11 February 2020), ImmuneSystemProcess-EBI-UniProt-GOA (11 February2020), MolecularFunction-EBI-UniProt-GOA (11 February 2020)), KEGG and Reactome. Fig. 1 | Histological analyses. a, Atlas of semi-quantitative synovial IHC scores for immune cells. b, Distribution of semiquantitative scores at baseline in all patients, individually shown in the y axis. The total on the x axis represents the sum of the individual scores (Immune score). c, Baseline semi-quantitative IHC scores, Krenn synovitis score ('Synovial score') and total Immune score in patients stratified according to 16 weeks CDAI50% response to rituximab (top) and tocilizumab (bottom). Two-sided Mann Whitney test. ns= p value >0.05. n = 161 patients. Boxplots showing median with first and third quartiles. Fig. 6 | Venn diagram showing overlap in genes between machine learning models and comparison with models built using only clinical and histological variables. a, Venn diagram showing the overlap in genes selected as features in optimal predictive models for prediction of rituximab and tocilizumab response at week 16 and refractory state (failure to respond to both rituximab and tocilizumab). b, Grid of plots showing the optimal predictive models for different treatment when using clinical and histological variables only. From top to bottom plots show: ROC curves for the best model on the test data (from outer-fold) set; ROC curves on the left-out (from inner-fold) set; and the variable importance when fit to the whole data set.