Introduction
The treatment pathway for rheumatoid arthritis (RA) includes conventional synthetic disease modifying anti-rheumatic drugs (csDMARD) [1]. Methotrexate (MTX) is the cornerstone csDMARD that is received by 60–88% of treated patients [2]. Despite its broad use, ∼30–40% of RA patients on MTX do not experience satisfactory disease control [3].
If symptoms are not effectively controlled by MTX, patients may receive another csDMARD (e.g. sulfasalazine) or may be escalated to receive an advanced targeted therapy (e.g. biologics) [4]. Inception cohort data have shown that early and effective treatment within the first 6 months of symptom onset can influence outcome decades later [5]. As response to first medication can have long-lasting effects and because the cost of advanced targeted therapies is reducing, being able to predict which patients would benefit most from MTX as a first-line therapy would be a significant clinical advancement.
Several clinical factors, including age, sex seropositivity and disability scores have shown moderate success in predicting MTX efficacy [6]; with female sex, seropositive disease and higher health assessment questionnaire scores (HAQ), a measure of functional disability, correlating with a less favourable response [7]. However, clinical models alone are not sufficiently predictive to guide clinical management, and they do not capture information on relevant biological pathways [7].
Genetic and genomic factors have been explored as potential predictors of MTX efficacy identifying biomarkers in transcriptomic [8], proteomic [9] and metabolomics modalities [10]. However, many previous studies have limited biomarker data assessment to a single timepoint, missing the opportunity to assess on-treatment change in biomarker profiles.
The current study aims to utilize comprehensive whole-blood RNA sequencing data collected at baseline and following 4 weeks on MTX and apply machine learning to identify pre- and early treatment biomarker profiles that are predictive of MTX treatment response, assessed at 6 months on therapy.
Methods
Patient cohort
The patients in this study were already taking part in the Rheumatoid Arthritis Medication Study (RAMS, Manchester NHS Research Ethics Committee reference 08/H1008/25), a national multicentre, longitudinal observational study in the UK that recruited patients with RA who were about to commence MTX therapy for the first time [7], and all patients provided written consent. RAMS patients included in the current study had a diagnosis of RA according to the American College of Rheumatology 1987 revised criteria for the classification of RA [11], were of European ancestry, were about to receive treatment with MTX for their RA symptoms and had available whole-blood samples for RNA extraction, and clinical data to address the study aims. In total, 100 MTX-treated patients were selected for this study for bulk transcriptomics at baseline and at 4 weeks follow-up.
Treatment response was assessed at 6 months after starting therapy using the European League Against Rheumatism (EULAR) criteria [12], where good responders were defined as patients with a change in disease activity score in 28 joints (DAS28-CRP) of >1.2 and a DAS28-CRP of <3.2 at 6 months, moderate responders were defined as patients with a change in DAS28-CRP of >0.6 and a DAS28-CRP of <1.2 at 6 months and non-responders were defined as patients with a change in DAS28-CRP of <0.6 at 6 months. Good and moderate responders were grouped together as ‘responders’ for analytical purposes.
RNA sequencing
RNA-Seq libraries were generated as previously described [13]. Briefly, libraries were enriched using oligo-dT beads from 500 ng of total RNA using the TruSeq stranded mRNA Library Prep Kit. Double-stranded cDNA was adenylated 3' and ligated to indexed adaptors with unique molecular identifiers (Integrated DNA Technologies) before being amplified by 15 PCR cycles.
Library quality was assessed using the Bioanalyzer DNA 7500 assay (Agilent Technologies, Santa Clara, CA, USA). Sequencing was performed on a NovaSeq 6000 (Illumina, San Diego, CA, USA) in paired-end mode with 2 × 151 bp read length. Image analysis, base calling and quality scoring were processed using RTA 3.4.4 analysis software, followed by the generation of FASTQ sequence files.
Reads were aligned to the GENCODE v39 hg38 human reference genome using STAR [14] and gene expression levels were quantified using RSEM [15].
Differential expression analysis
Genes with cumulative counts under 500 across samples and those expressed in <20% of samples were excluded. The glmmSeq R package [16] was used for differential gene expression analysis using a generalized linear mixed model with negative binomial distribution, including patient variability as a random effect. DESeq2 [17] was used to estimate dispersion parameters.
Statistical models were adjusted for sex, baseline DAS28-CRP score and an interaction term was specified between sampling timepoint and treatment response. Sex was chosen as a covariate due to the female predominance of the disease being treated, and DAS28-CRP was included because preliminary cohort comparisons indicated a significant difference in baseline disease activity between responders and non-responders. The resulting model is:
Where yij is the gene expression for patient i at observation j; β0 is the intercept; β1 and β2 are fixed-effect coefficients for the covariates sex and baseline DAS28-CRP scores, respectively; β3, β4, and β5 are fixed-effect coefficients for response, timepoint and their interaction, respectively; and ui∼N(0,σu2) is the random effect term for patient variability. Sex and baseline DAS28-CRP were included as covariates in the models to account for potential confounding. Baseline DAS28-CRP differed between responders and non-responders, while sex was included a priori as a clinically relevant variable in RA.
A false discovery rate (FDR) below 0.05 and a fold-change exceeding 1.25 in either the interaction or treatment response identified genes as differentially expressed.
Machine learning
DESeq2’s variance stabilizing transformation (VST) was applied for normalizing variance across genes prior to multivariate analyses. Features going into the model included normalized gene expression values, sex and baseline DAS28-CRP scores. Training of Gradient Boosted Trees Classifier models was carried out using the XGBoost library [18].
XGBoost was chosen due to its ability to handle class-imbalanced dataset with weights on classes, high suitability for tabular data [19] and compatibility with SHapley Additive exPlanations (SHAP) [20], which is a method to interpret machine learning model. Nested stratified 10-fold cross-validation was used to avoid overfitting of model performance. To reduce the high dimensionality of gene expression data, glmmSeq was used to select differentially expressed genes according to the criteria described above, using the repeated-measures model described in the preceding section rather than separate timepoint-specific analyses.
Gene selection was performed within classifier training to prevent information leakage between training and testing datasets and to minimize overfitting. Feature (gene) selection and model performance were evaluated using nested cross-validation, with an outer 10-fold cross-validation loop for model evaluation and an inner 5-fold loop for hyperparameter optimization. Within each outer training fold, feature selection was performed independently, allowing the stability of selected genes to be assessed across resampled training subsets while preserving independence of the held-out test data.
Hyperparameter tuning of the XGBoost model was carried out using scikit-optimize’s Bayes search CV [21] within the inner cross-validation loop using a search over predefined parameter ranges. The parameters tuned included the number of trees (n_estimators), maximum tree depth (max_depth), learning rate, row subsampling ratio (subsample), column subsampling ratio per tree (colsample_bytree), minimum loss reduction for node splitting (gamma) and L1 and L2 regularization parameters (reg_alpha and reg_lambda).
The search ranges were: n_estimators (100–2,000), max_depth (1–10), learning_rate (0.001–1.0), subsample (0.05–1.0), colsample_bytree (0.05–1.0), gamma (0–0.7), reg_alpha (0–1) and reg_lambda (0–1). The optimal parameter combination was selected based on the mean area under the receiver operating characteristic curve (ROC–AUC) across the inner cross-validation folds. Three models were explored: (i) baseline only data; (ii) 4-week only; and (iii) inputs from both time-points (baseline and 4-week).
To benchmark predictive performance against clinical characteristic data, an additional clinical characteristics only XGBoost classifier was trained, with hyperparameters optimized using 10-fold cross validation. Following fitting of the model, the model was interpreted using SHAP values averaged over 1,000 initializations of gradient boosted trees to account for variability in this stochastic process. Finally, features with mean absolute SHAP values of at least 0.05 were considered the most important features and were used to train the final model, and the same set of features were used for functional enrichment analysis using Metascape [22].
Cell deconvolution analysis
To assess whether prediction is driven by immune-cell composition rather than gene regulation within cells, cell-type proportions were estimated from RNA-seq data using CIBERSORTx [23] with the LM22 reference panel. Downstream transformation, statistical testing and multiple-testing correction were performed according to our previously published workflow [13].
Results
Following data quality control 16 548 genes were available for analysis in 100 patients. Patients were either categorized as responders (n = 69) or non-responders (n = 31) following 6 months of MTX therapy. The baseline characteristics of the patients are presented in Table 1.
Baseline characteristics were broadly comparable between responders and non-responders, with no statistically significant differences observed for the variables listed in Table 1, except for baseline DAS28-CRP.
Patient baseline characteristics: differential expression
Ninety-three genes were differentially expressed between responders and non-responders across both timepoints. Eighty-one of these were not statistically significant in the interaction term; that is, they were not dependent on time. Arachidonate 15-Lipoxygenase, Type B (ALOX15B) was the most upregulated gene in responders in the mixed effect model, which accounted for both timepoints, showing a fold change of 1.7.
This was followed by Oleoyl-ACP Hydrolase (OLAH) and CUB and Sushi Multiple Domains 1 (CSMD1) with fold changes of 1.56 and 1.58, respectively (Table 2). Nebulette (NEBL), associated with actin cytoskeletal organization, was the most downregulated gene in responders, with a fold change of 0.52 (Table 2). Twelve differentially expressed genes were dependent on time (Figure S1).
Expression of OLAH, IGHV3-72 and AMPH was reduced in responders by 0.76, 0.76 and 0.68-fold over the 4 weeks, respectively (Table 2, Fig. 1). Enrichment analysis with Metascape revealed that the majority of response associated genes were enriched in adaptive immune response and fatty acid biosynthetic process (Figures S2–S4).
Machine learning classifier of methotrexate response
A machine learning classifier for predicting treatment response outcome built using only clinical characteristics was first explored (not shown). Firstly, a basic clinical model including patient sex and baseline DAS28-CRP was developed and achieved a mean internal cross-validation AUC of 0.68 (± 0.13). A second, more complex model including all clinical information listed in Table 1 achieved similar performance with a mean internal cross-validation AUC of 0.72 (± 0.11).
Transcriptomics model development retained patient sex and baseline DAS28-CRP to maintain a parsimonious model and reduce the risk of overfitting. The machine learning model that included gene expression and clinical information incorporating both timepoints performed less well than single timepoint models, with an AUC of ∼0.84. We therefore prioritized the models built with gene expression and clinical information that require less sampling.
These models, stratified by individual timepoints, achieved an AUC of 0.89 (± 0.09) and 0.90 (± 0.08) for the baseline and 4-week classifier, respectively (Fig. 2). Both models also achieved good performance in other classification metrics (Table 3).
Using mean absolute SHAP values of 0.05 as cut-off, the baseline model had 33 features surpassing this threshold and the 4-weeks model had 54. One thousand initialization of gradient boosted trees were evaluated for SHAP values, and their standard deviations were computed and shown in Fig. 3; the standard deviations are larger than the mean for some features, particularly seen in the 4-week model.
Among the important set of genes, six genes at baseline and 15 genes at week 4 have unknown functions, with annotations available only through their Ensembl IDs or as known long non-coding RNAs.
At baseline, the three most informative genes with known annotations are Lipocalin-12 (LCN12), Caveolin-1 (CAV1) and Galactosidase Beta 1 Like (GLB1L). The first two genes demonstrated significant downregulation in responders, with CAV1 downregulated by 0.92, LCN12 by 0.79, but GLB1L was upregulated by 1.26. CAV1 is a scaffolding protein that plays a critical role in membrane structure and signal transduction.
It has been implicated in cancer, cardiovascular and pulmonary diseases due to its tumour suppressor functions and regulatory influence on cellular pathways [33]. LCN12 belongs to the lipocalin family and is involved in the transport of hydrophobic molecules such as steroids and has been associated with inflammation and cancer [34]. GLB1L encodes a predicted β-galactosidase-like protein within the glycosyl hydrolase family, suggesting potential involvement in glycan processing pathways, although its precise biological role remains poorly characterized [35].
Other notable genes in the baseline model include: lipid transport protein prostaglandin D synthase (PTGDS), a potential prognostic marker in cancer that is associated with activation of immune cells [36]; TNF superfamily member 15 (TNFSF15), a member of tumour necrosis factor family that is inducible by TNF, is expressed on several immune cells and is abnormally expressed in inflammatory disease including RA [37]; and HERV-H LTR-Associating 2 (HHLA2) an immune checkpoint molecule that has immune inhibitory and stimulatory functions, depending on the type of receptor it binds to [38].
At week 4, the three most informative genes with known annotations were Scavenger Receptor Family Member Expressed on T Cells 1 (SCART1), sulfotransferase family 1A member 2 (SULT1A2) and RAR-Related Orphan Receptor C (RORC). SCART1 and RORC exhibited a 21% downregulation in responders, while SULT1A2 is upregulated by 32%. SCART1 is a gene that likely functions as a co-receptor on αβ and γδ T cells, playing a role in immune system regulation; however, its specific functions and disease associations remain poorly characterized.
SULT1A2 encodes a cytosolic sulfotransferase belonging to the SULT1 family, which catalyses sulphf conjugation reactions in phase II metabolism; however, the precise biological role of SULT1A2 remains less well defined than that of better studied sulfotransferases [39]. RORC is a nuclear receptor that functions as a DNA-binding transcription factor essential for lymphoid organogenesis and thymopoiesis. Dysregulation of RORC has been associated with immunodeficiency and inflammatory diseases [40].
Both gene sets were statistically significantly enriched in immune system processes.
Estimated immune-cell proportions derived from CIBERSORTx were not significantly different between responders and non-responders, with all comparisons yielding unadjusted P-values >0.05.
Discussion
In this study, we performed whole-blood RNA-sequencing at pre-treatment and following 4-weeks on treatment to develop a classifier of MTX response by 6 months and identify gene biomarkers in responders and non-responder patient groups, where responders are patients with EULAR treatment response classification of good and moderate. This grouping was selected to capture informative biomarker profiles prior to any improvement on treatment. Overall, genes found to be associated with response either through statistical modelling or machine learning model were enriched in adaptive immune responses.
The machine learning models in this study achieved a higher AUC of ∼0.89 when individual timepoints were analysed, an improvement over previous blood-based studies using microarrays [8] or in smaller sample numbers [41].
Among the differentially expressed genes, several have established associations with RA. For instance, ALOX15B is involved in the metabolism of polyunsaturated fatty acids and has been implicated in inflammatory processes and is upregulated in RA patients [42]. OLAH and Versican (VCAN) are involved in fatty acid biosynthesis and cell adhesion in RA patients within bone marrow and synovial fibroblasts, respectively [25].
All three genes were upregulated in responders with expressions decreasing in the post treatment sample (Supplementary Fig. S5). Other genes that were upregulated in responders but decreased following treatment with methotrexate included IGHV3-72 and AMPH.
The machine learning approach revealed that the baseline model was the most informative and stable, based on dispersion of SHAP values. The model trained on the 4-week timepoint resulted in a lower score on all the metrics as well as high standard deviation on the SHAP values, indicating potential overfitting and less generalizable results. The difference in genes identified between the statistical model and machine learning model were observed with high effect sizes or low P-value not translating to high mean absolute SHAP values (Supplementary Table S1).
Despite this, the genes identified by both methods are implicated in various processes with relevance to RA including inflammation and atherogenesis: ALOX15B [43]; complement activation: CSMD1 [44]; cell proliferation, migration and adhesion: SDC1 [31] and proliferation and cytokine production in T cells: HHLA2. HHLA2 can either bind to the receptor KIR3DL3 or TMIGD2. TMIGD2 is expressed on various subsets of T cells and Natural Killer (NK) cells and interaction with HHLA2 [38] results in enhanced proliferation and cytokine production.
Conversely, binding of HHLA2 to KIR3DL3 results in inhibition of T cells and NK cell pathways [38]. Given the dual role of HHLA2, the importance of checkpoint molecules in both cancer and inflammatory diseases [45], and the identification of HHLA2 in both analyses presented here, this molecule deserves further investigation in independent experiments.
The identified genes also revealed a theme of apoptosis. For example, Klotho (KL) and VCAN were up-regulated in responders with 1.27-fold and 1.4-fold change, respectively, whilst NSG1 was down-regulated in responders with 0.76-fold change. They were ranked 26, 28 and 9 based on baseline classifier’s mean absolute SHAP values.
NSG1 is a known mediator of p53 inducing apoptosis [46], while KL has been reported to promote apoptosis in cancer cells [47]; in contrast, VCAN has been shown to exert anti-apoptotic effects, particularly under oxidative stress conditions [48]. The involvement of these genes may indicate a possible link between MTX use and apoptotic processes, aligning with the drug’s known mechanisms of action [49]. However, experiments in enriched cell types that are relevant to RA are required to investigate this further.
The current study has several strengths and limitations that should be acknowledged. The RAMS study is observational and is vulnerable to biases and confounding that may only partially be controlled for during data analysis. Common confounding would include patients having different baseline DAS28, which were accounted for during analysis.
Further, the use DAS28-CRP composite score is likely to capture information beyond inflammation or mechanisms that have a clear biological basis. Although CDAI would provide an alternative outcome measure that does not incorporate CRP, its use was not possible because of the high proportion of missing data for the physician global assessment component. The feature selection strategy employed may not have captured the full spectrum of gene interaction because the selection was carried using univariate analysis prior to machine learning modelling.
While this approach reduced the dimensionality of gene-expression matrix and helped limit overfitting in the context of a modest sample size, it may have excluded genes with predictive value only in combination with other genes. Future studies with larger cohorts could compare alternative feature-selection strategies, including variance-based filtering, penalized regression and embedded feature selection methods. Strengths include the comprehensive RNA sequencing approach across multiple time points, allowing the change in expression profile following treatment to be observed.
The machine learning pipeline employed robust internal validation coupled to the powerful gradient boosting tree algorithm that can incorporate weights into the model to account for class imbalance. Although we mitigated overfitting with the modest sample size by cross-validation, the lack of validation in an independent external cohort is a limitation that prevents assessment of the generalizability of the model.
Conclusion
In conclusion, this study provides valuable insights into the genomics underpinnings of MTX response in RA, offering a basis for future research aimed at improving personalized treatment strategies by including blood-based biomarker profiling. Despite the limitations outlined above, the study’s novel approach and significant findings contribute to the growing body of knowledge in RA treatment and precision medicine, but findings require independent replication.
Supplementary Material