Multiomics and Machine Learning Accurately Predict Clinical Response to Adalimumab and Etanercept Therapy in Patients With Rheumatoid Arthritis

Abstract

Objective: To predict response to anti-tumor necrosis factor (anti-TNF) prior to treatment in patients with rheumatoid arthritis (RA), and to comprehensively understand the mechanism of how different RA patients respond differently to anti-TNF treatment.

Methods: Gene expression and/or DNA methylation profiling on peripheral blood mononuclear cells (PBMCs), monocytes, and CD4+ T cells obtained from 80 RA patients before they began either adalimumab (ADA) or etanercept (ETN) therapy was studied. After 6 months, treatment response was evaluated according to the European League Against Rheumatism criteria for disease response. Differential expression and methylation analyses were performed to identify the response-associated transcription and epigenetic signatures. Using these signatures, machine learning models were built by random forest algorithm to predict response prior to anti-TNF treatment, and were further validated by a follow-up study.

Results: Transcription signatures in ADA and ETN responders were divergent in PBMCs, and this phenomenon was reproduced in monocytes and CD4+ T cells. The genes up-regulated in CD4+ T cells from ADA responders were enriched in the TNF signaling pathway, while very few pathways were differential in monocytes. Differentially methylated positions (DMPs) were strongly hypermethylated in responders to ETN but not to ADA. The machine learning models for the prediction of response to ADA and ETN using differential genes reached an overall accuracy of 85.9% and 79%, respectively. The models using DMPs reached an overall accuracy of 84.7% and 88% for ADA and ETN, respectively. A follow-up study validated the high performance of these models.

Conclusion: Our findings indicate that machine learning models based on molecular signatures accurately predict response before ADA and ETN treatment, paving the path toward personalized anti-TNF treatment.

Figure 1

Flow chart showing the study methods. A, Blood samples were obtained from patients with rheumatoid arthritis (RA) at baseline, and patients were treated with subcutaneous adalimumab (ADA) or etanercept (ETN) for 6 months. Peripheral blood mononuclear cells (PBMCs) were isolated for RNA sequencing and DNA methylation profiling. CD4+ T cells and CD14+ monocytes were then isolated for RNA sequencing. Patients were classified as responders or nonresponders according to European League Against Rheumatism criteria at the end of month 6. Patients who did not respond to ADA were switched to ETN, and patients who did not respond to ETN were switched to ADA, for 6 months and treatment responses were observed. B, RNA sequencing data sets and DNA methylation data sets were used to build machine learning models to predict treatment responses at baseline.

Figure 2

Differences between gene expression in PBMCs from responders (R) versus nonresponders (NR) to ADA or ETN. A and B, Heatmaps showing the top 100 differentially expressed genes (DEGs) in PBMCs from responders versus nonresponders to ADA (A) and ETN (B). C and D, Expression of example DEGs in responders and nonresponders to ADA (C) and ETN (D). Data are shown as box plots in the style of Tukey. Each box represents the 25th to 75th percentiles. Lines inside the boxes represent the median. Whiskers outside the boxes represent the lowest and highest data points excluding the outliers. The outliers are determined by 1.5‐interquartile range criteria. Circles represent individual patients. E and F, Multidimensional scaling plots based on gene expression data showing clusters of patient response to ADA (E) and ETN (F). G, Venn diagram showing the overlap in DEGs between patients receiving ADA and patients receiving ETN. H, Scatterplot showing differences in fold change (FC) in RNA expression for DEGs in the ADA cohort and the ETN cohort. * = P < 0.05; ** = P < 0.01. CPM = count per million; NS = not significant (see Figure 1 for other definitions). Color figure can be viewed in the online issue, which is available at http://onlinelibrary.wiley.com/doi/10.1002/art.41516/abstract.

Figure 3

Differentially methylated positions (DMPs) associated with response to ADA and ETN in PBMCs. A and B, Manhattan plots showing the −log10‐transformed P values for CpG sites in the DNA methylation profiling for response to ADA (A) and ETN (B). Gene names are shown for genes in which any CpG site reached a significance level of P < 10⁻⁴. The corresponding probes are shown in Supplementary Figure 2, available on the Arthritis & Rheumatology website at http://onlinelibrary.wiley.com/doi/10.1002/art.41516/abstract. C and D, Multidimensional scaling plots, based on DNA methylation data, showing clusters of patients classified as responders or nonresponders to ADA (C) and ETN (D). E and G, Distribution of gene region features where DMPs associated with response to ADA (E) and to ETN (G) were located. TSS1500 = 200–1,500 bases upstream of the transcription start site (TSS); TSS200 = 0–200 bases upstream of the TSS; 5′UTR = 5′‐untranslated region; ExonBnd = exon boundary; IGR = intergenic region. F and H, Distribution of area related to CpG island where DMPs associated with response to ADA (F) and ETN (H) were located. island = CpG island; shore = 0–2 kb from island; shelf = 2–4 kb from island; opensea = the rest of the area. See Figure 1 for other definitions. Color figure can be viewed in the online issue, which is available at http://onlinelibrary.wiley.com/doi/10.1002/art.41516/abstract.

Figure 4

Differences between gene expression in monocytes and CD4+ T cells from responders (R) versus nonresponders (NR) to ADA and ETN. A and B, Heatmaps showing the top 100 differentially expressed genes (DEGs) in monocytes (A) and CD4+ T cells (B) from responders versus nonresponders to ADA. C, KEGG pathways enriched by the DEGs associated with ADA response in monocytes and CD4+ T cells. D, Expression of ADA response–associated DEGs involved in RA and the tumor necrosis factor (TNF) signaling pathway in CD4+ T cells. E and F, Heatmaps showing the top 100 DEGs associated with ETN response in monocytes (E) and CD4+ T cells (F). G, KEGG pathways enriched by the DEGs associated with ETN response in monocytes and CD4+ T cells. H, Expression of ETN response–associated DEGs involved in FoxO signaling, NOD‐like receptor signaling, and JAK/STAT signaling pathways in CD4+ T cells. I, Venn diagrams showing the overlap between ADA response–associated DEGs and ETN response–associated DEGs in monocytes and CD4+ T cells. J and K, Venn diagrams showing the overlap between ADA response–associated DEGs (J) and ETN response–associated DEGs (K) in monocytes, CD4+ T cells, and PBMCs. In D and H, data are shown as box plots in the style of Tukey. Each box represents the 25th to 75th percentiles. Lines inside the boxes represent the median. Whiskers outside the boxes represent the lowest and highest data points excluding the outliers. The outliers are determined by 1.5‐interquartile range criteria. Circles represent individual patients. * = P < 0.05; ** = P < 0.01; *** = P < 0.001. CPM = count per million (see Figure 1 for other definitions). Color figure can be viewed in the online issue, which is available at http://onlinelibrary.wiley.com/doi/10.1002/art.41516/abstract.

Figure 5

Machine learning models of the prediction of RA patients’ response to ADA and ETN therapy. A and B, Accuracy, true‐negative rate (TNR), and true‐positive rate (TPR) of machine learning models based on gene expression signatures (CD14, CD4, and PBMC RNA) and DNA methylation signatures (PBMC DNA) for the prediction of 6‐month response to ADA (A) and ETN (B). Bars show the mean ± SEM from cross‐validation analysis. C and D, Machine learning prediction of 6‐month response to ADA and ETN treatment using the best model based on gene expression data (C) and DNA methylation data (D). R = responder; NR = nonresponder; AR = responder to ADA only; ER = responder to ETN only; DR = double responder (response to both ADA and ETN). E, Validation of machine learning prediction by 6‐month‐follow‐up drug‐switched study (see Figure 1). The last 2 columns are predictions made by the best machine learning models based on gene expression data and DNA methylation data, respectively. Correct predictions are shown in red. See Figure 1 for other definitions. Color figure can be viewed in the online issue, which is available at http://onlinelibrary.wiley.com/doi/10.1002/art.41516/abstract.