Identification of Three Rheumatoid Arthritis Disease Subtypes by Machine Learning Integration of Synovial Histologic Features and RNA Sequencing Data

Patient data.

An overview of the study design is presented in Figure 1. We enrolled 123 consecutive patients with RA undergoing arthroplasty at the Hospital for Special Surgery (HSS) in New York who had met either the American College of Rheumatology (ACR)/European League Against Rheumatism 2010 classification criteria for RA (18) and/or the ACR 1987 criteria for RA (19). In addition, synovial samples from 6 patients with OA were included in the analyses.

The serum status of ACPAs (measured with the anti–cyclic citrullinated peptide antibody [anti-CCP] test) and RF was coded as negative, low or medium positive (1–3 times the upper limit of normal), or high positive (>3 times the upper limit of normal). Pain was assessed by asking patients the question, “How much pain have you had because of your condition over the past week? Please indicate how severe your pain has been on a scale of 0–10.” This study was approved by the HSS Institutional Review Board (approval no. 2014–233), the Rockefeller University Institutional Review Board (approval no. DOR0822), and the Biomedical Research Alliance of New York (approval no. 15–08-114–385). All participating patients provided their signed informed consent.

Sample processing.

Adjacent areas of synovial tissue were placed into OCT frozen blocks for histology and processed for cell dissociation. For dissociation, the tissue was dissected and treated with Liberase TL (100 μg/ml; Roche) and DNase I (100 μg/ml; Roche) at 37°C for 15 minutes. Red blood cells were lysed and the remaining cells were washed and filtered through 70-μM strainers. Synovial cells were lysed in RLT (Qiagen) with 1% β-mercaptoethanol, and then stored at −80°C.

RNA sequencing.

RNA was extracted using a Qiagen RNeasy Mini kit (no. 74104). Libraries were prepared using Tru-seq messenger RNA (mRNA) Stranded Library kits; 50-bp, paired-end reads were sequenced on a HiSeq2500 platform. Reads were aligned to hg19 using STAR (20). Samples with >0.1% globin mRNA were excluded from further analysis to prevent confusion between infiltrating and contaminating hematopoietic cells. Forty-five samples of sufficient quality were processed in 3 separate batches. To account for any batch effects, we used ComBaT in the Bioconductor SVA package (21) and DESeq2 (22) to normalize the data.

Clustering of RNA-seq data.

RNA-seq data were clustered using consensus clustering, an iterative clustering method in which small data perturbations (using the rnorm function in the basic stats package in R) are introduced at each clustering iteration. The clustering iterations offer the opportunity to derive a statistically robust data partitioning, as opposed to a one-time, single clustering computation. We used the k-means clustering method, which typically does not converge to the same solution when clustering the same data multiple times (contrary to agglomerative approaches, such as hierarchical clustering). We clustered the 45 RNA-seq samples 1,000 times into k = 2, k = 3, k = 4, and k = 5 clusters. The clustering iterations were then used to generate likelihood scores representing the frequency of a pair of samples clustering together (number of times the pair of samples cluster together divided by the number of clustering iterations). This score ranges from 0 to 1, with scores of 1 indicating co-clustering of the 2 samples.

Differential gene expression analysis.

We used the DESeq2 Bioconductor package (22) to identify differentially expressed genes (DEGs) across the synovial subtypes identified via clustering, at an adjusted P (P_adj) value of 0.01. We eliminated chromosome X/Y genes to remove sex biases, as well as IgG variable (V D, and J) genes, since the individual V-, D-, and J-gene segments are counted as individual immunoglobulin genes, thereby overrepresenting their abundance. Immunoglobulin constant regions were retained to ensure representation of immunoglobulin gene expression.

Pathway analysis.

DAVID analytic modules (version 6.8) (23) with default parameters were used to identify functional annotation clusters of genes. Genes found to be differentially expressed in the high inflammatory subgroup compared to the other subgroups and low inflammatory subgroup compared to the other subgroups, with a P_adj value of less than 0.05 and a log₂-transformed fold change of greater than 0.5, were used to identify enriched pathways. Enrichment scores were the modified P values from Fisher’s exact test based on data from the gene-enrichment analysis.

Inference of cell types in gene expression clusters.

In order to deconvolute the cellular composition of the 3 subtypes in our data, we used an algorithm called CIBER-SORT (Cell-type Identification by Estimating Relative Subsets of Known RNA Transcripts) (24,25), a machine learning system trained on the profiles of 22 distinct leukocyte data sets over 547 genes, with defined “barcodes” of gene expression signatures to distinguish cell types.

Predicting RNA-seq subtypes from histologic features.

We used histologic scores as modeling features in a standard SVM (Bioconductor package e1071) to predict the 3 RNA-seq subtypes, using leave-one-out cross validation to train and predict transcriptomic subtypes. Weights attributed to each feature were extracted from an SVM. All histologic features were presented as binary vectors—histologic features having more than 2 categories were converted to a binary representation, so that all features were equally represented in the model.

Histologic scoring.

In total, 129 synovial samples were preferentially obtained from grossly inflamed (dull and opaque) synovium. If no inflammation was apparent, samples were obtained from standard locations: the femoral aspects of the medial and lateral gutters, and the central supratrochlear region in the suprapatellar pouch. Tissue samples were snap-frozen and stained with Harris’ modified hematoxylin solution and eosin Y (both from Sigma-Aldrich). Twenty histologic features (synovial lining hyperplasia, lymphocytes, plasma cells, Russell bodies, binucleated plasma cells, fibrin, synovial lining giant cells, sublining giant cells, fibrosis, detritus, necrosis, granulation tissue, neutrophils, mast cells, eosinophils, synovial chondrometaplasia, germinal centers, mucin, infarction, and vascularity) were scored using a systematic approach (as outlined in detail in Supplementary Methods, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract and also available at www.hss.edu/pathology-synovitis). In total, 129 synovial samples were scored by a single pathologist (EFD), and 40 samples were scored again by a second pathologist (EG) to determine interrater reliability.

ACPA array.

As previously described (24), serum levels of antibodies targeting 38 putative RA-associated autoantigens were measured using a custom bead-based immunoassay. Serum samples were diluted at a 1:30 ratio in a proprietary dilution buffer (Bio-Rad), mixed with the antigens conjugated to spectrally distinct fluorescent microspheres (Bio-Rad), and then incubated with an anti-human IgG antibody conjugated to phycoerythrin (Jackson ImmunoResearch). The resulting fluorescence intensities were quantitated using a Luminex 200 System.

Statistical analysis.

Kappa statistics were used to assess the interrater reliability of the histologic scores sourced from the 2 pathologists. For variables with ≥3 response levels, weighted kappa statistics were calculated to give credit for differences that were close, as opposed to treating any difference as the same. Interrater reliability was scored as follows: none to slight, κ = 0.01–0.2; fair, κ = 0.21–0.40; moderate, κ = 0.41–0.60; substantial, κ = 0.61–0.80; and almost perfect, κ = 0.81–1.00 (25). The Jonckheere-Terpstra 2-sided test for trend was used to compare pathway enrichment scores from the gene set variation analysis data across the 3 synovial subtypes. Kruskal-Wallis test with Dunn’s test for multiple comparisons was used to compare CIBERSORT RNA scores and clinical features among the 3 synovial subtypes. The chi-square test was used to detect differences in binary data above those expected by chance among the synovial subtypes. Analysis of variance with Tukey’s test for multiple comparisons was used to compare the log₂-transformed mean fluorescence intensities of ACPA expression among the 3 synovial subtypes. Spearman’s correlation coefficients were calculated to assess correlations between the severity of pain and gene expression levels across the 3 synovial subtypes.

Clinical characteristics.

The clinical characteristics of the 123 patients with RA are presented in Supplementary Table 1 (available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract). The majority of the RA patients were female. Among the 123 patients, 47% were seropositive for RF and 50% were seropositive for anti-CCP. Although the mean disease duration of the RA patients was 14 years, the mean Disease Activity Score in 28 joints was 3.8 (indicative of moderate disease activity), and 53% of the patients were treated with a biologic agent just prior to surgery. The affected joints included the hips in 40% of patients, knees in 57% of patients, and shoulders in 3% of patients.

Histologic features of the extracted synovium.

We first sought to determine the prevalence and feasibility of scoring the candidate histologic features of the RA and OA synovium. Fourteen of the 20 features were observed in >5% of synovial samples (Figure 2A).

We next evaluated the interrater reliability for the scoring of these histologic features on a subset of samples (Figure 2B). Representative images of the features with a frequency of >5% in the synovial samples and having at least fair interrater reliability are presented in Figures 2C and D. Features that were seen in at least 5% of samples and that had at least fair interpathologist reliability included plasma cells, binucleated plasma cells, Russell bodies, giant cells, neutrophils, fibrin, mucin, detritus, synovial lining hyperplasia, and lymphocytes.

Three distinct subtypes of synovium identified by RNA-seq data analysis.

Independent of the histologic analysis, consensus clustering of the top 500 most variable genes expressed in the 45 synovial samples (39 from patients with RA and 6 from patients with OA) identified 3 gene expression subtypes. Figure 3A depicts the pairwise likelihood scores for the likelihood of samples clustering together, with the 45 samples being constrained to k = 2, k = 3, and k = 4 clusters. The red boxes in Figure 3A represent the co-clustering samples (likelihood scores of 1), blue represents the samples that never coclustered (scores of 0), and white represents inconsistent co-clustering patterns (scores of 0.5). When the samples were partitioned into 2 or 3 clusters, the likelihood score was crisp and well defined, with the sharp red and blue colors suggesting that the majority of the likelihood scores approached 1 or 0. However, when the samples were clustered into 4 or more groups, significantly less consistency was observed. This visualization is a statistically robust confirmation that the 45 RNA-seq samples formed, at most, 3 distinct subgroups.

We also explored this analysis using differing numbers of variable genes (see Supplementary Figure 1, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract) and identified very minor shifts across the 3 consensus clustering subtypes, lending confidence to the k = 3 sample partitioning. Principal components analysis of the samples using the top 500 most variable genes demonstrated agreement with the findings from the RNA-seq analyses that revealed the 3 synovial subtypes (Figure 3B), thus validating the clustering analysis.

Differential gene expression patterns between the 3 synovial subtypes.

We identified 6,582 transcripts as DEGs (P_adj < 0.01) across the 3 synovial clusters. Comparisons of each cluster to the others are presented in Supplementary Tables 2, 3, and 4 (available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract). More than one-half of these DEGs separated the low inflammatory from the high inflammatory subtype, while the majority of the remaining genes separated the mixed from the high inflammatory subtype. The mixed subtype shares features with both the high and low inflammatory subtypes. Since it is a blend of the 2 subtypes, there are very few genes that are unique to the mixed group. Heatmaps of the top 1,000 and top 50 most variable genes (as ranked by DESeq2 analysis based on the P_adj value output) depict the main gene expression patterns segregating the 3 subtypes (Figure 3C). The largest gene block displayed an increasing expression pattern as the categories progressed from low to mixed to high inflammatory subtype (see Supplementary Table 5, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract).

Using DAVID gene ontology analysis, we performed functional annotation clustering on the genes with increased expression in the high inflammatory subtype compared to the others, and found that the highest enrichment scores corresponded to the pathways of immunity, immune cell signaling (such as SH2 and SH3 domains and kinases), and immunoglobulins, as well as chemokines and cytokines (Figure 4A) (see also Supplementary Table 6, http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract).

A smaller block of genes showed higher expression in the low inflammatory and mixed subtypes than in the high inflammatory subtype (Figure 3C) (see also Supplementary Table 7, http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract). Functional annotation clustering of genes with increased expression in the low inflammatory subtype compared to the other subtypes identified enrichment in pathways that included glycoproteins (Figure 4B) (see also Supplementary Table 8, http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract). These genes included markers of lining layer fibroblasts such as CD55, transforming growth factor β (TGFβ) superfamily genes (TGFBR1, TGFBR2, TGFB2, and BMP6), and genes involved in the regulation of extracellular glycoproteins such as GCNT3 (a glycosyltransferase that mediates core O-glycan branching, a critical step in mucin synthesis) (26–30). Another enriched pathway in the low inflammatory subtype included cell adhesion (Figure 4B). This pathway includes protocadherin gamma genes, a family of cell-adhesion proteins predominantly expressed in neuronal synapses that are critical for axon targeting and neuron survival (31–34).

Other significantly enriched pathways were those involved in neuronal pain processing. For example, ADRA1B (α−1 adrenergic receptor) is up-regulated on nociceptive nerve fibers that survive peripheral injury (35–37). We also analyzed expression of genes previously identified to discriminate lymphoid, myeloid, fibroid, and low inflammatory synovial subtypes (5). Samples from our high inflammatory subtype demonstrated elevated expression of myeloid genes as well as lymphoid signature genes (see Supplementary Figure 2, http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract). The previously described “fibroid” synovial subtype was enriched in TGFβ signaling pathway genes, with a relative paucity of inflammatory gene expression (5). We confirmed that these genes were also overexpressed in our low inflammatory subtype (Supplementary Figure 2).

CIBERSORT deconvolution of synovial gene expression subgroups.

Identifying the cellular source of gene expression variation can be challenging in samples containing various cell types, due to differences both in proportions and in total cell quantities. To better characterize the cell types responsible for gene expression differences among the synovial tissue samples, we applied the machine learning framework of CIBERSORT to estimate cell type composition. We compared inferred leukocyte frequencies using CIBERSORT analysis of gene expression profiles across the 3 synovial subtypes. The clusters were characterized by a successive increase in the proportion of hematopoietic cells between the low (23%), mixed (32%) and high inflammatory (68%) subgroups (Figure 4C). The high inflammatory subgroup harbored significantly increased inferred fractions of neutrophils as well as M2, M1, and M0 macrophages, monocytes, T follicular helper cells, memory activated and resting CD4 T cells, CD8 T cells, and plasma cells.

Predicting gene expression subtypes from histologic features using machine learning.

One of our main goals was to determine the synchrony between synovial histologic features and their genomic subtypes, the existence of which enables a cheaper histology-based approach to characterizing synovial tissue. To this end, we implemented a leave-one-out cross-validation SVM classification system to predict synovial genomic subtypes for our 45 samples, using a binary representation for the histologic scores as training features and genomic subtypes as training labels for the model. The predictive power of this system was evaluated by measuring the area under the receiver operating characteristic curve (AUC) (Figure 5A). Models separating the high inflammatory subtype from the others and those separating the low inflammatory subtype from the others performed best (AUCs of 0.88 and 0.71, respectively), whereas models separating the mixed subtype from the others were harder to predict (AUC of 0.59).

We compared the predicted synovial subtypes to those assigned by RNA-seq clustering, and found agreement with the RNA-seq data in 39 of the 45 samples. Histologic feature weights were extracted from the modeling system (as shown in Figure 5B), and we observed that, in general, the most discriminatory features also had stronger interrater reliability (Figure 2B).

After observing satisfactory training and testing results across the 45 RNA-seq samples, we computed 2 SVMs using all 45 samples: one SVM was used to discriminate the low inflammatory subtype from the others, and one was used to discriminate the high inflammatory subtype from the others. We used these models to predict genomic subtypes of the remaining 82 synovial samples for which we had complete histologic and clinical laboratory data available, but not gene expression data. We compared the distributions of histologic features across the RNA-seq samples from 39 RA patients classified by gene expression clustering to the distributions of histologic features across the 82 samples classified by the histologic scoring algorithm. The distributions of histologic features were similar whether the samples were classified by gene expression clusters or by the histologic feature algorithm (Figure 5C). For example, lymphocyte scores were always highest in the high inflammatory subtype, regardless of how the subtypes were classified. This is an expected result and serves to validate our approach.

We also compared the distributions of clinical laboratory features between the 39 RA patients classified by gene expression clusters and the 82 RA patients classified by the histologic scoring algorithm (Figure 5D), and again found consistent patterns. Given that clinical features were not part of the development of the algorithm, this consistent distribution of clinical features demonstrates that both classification methods showed associations in meaningful ways with clinical laboratory features. For example, low ESR values were most commonly present in the low inflammatory subtype and rarely present in the high inflammatory subtype, whether classified by gene expression or the histologic scoring algorithm. This consistent partitioning of an independent set of clinical features provides additional confirmation for the histologic scoring algorithm.

Association of synovial subtypes with clinical features.

We next compared the clinical features of all 123 RA patients who were classified using the histologic scoring algorithm. First, we compared the histologic features across the 3 subtypes and found that the levels of lymphocytes, plasma cells, synovial lining hyperplasia, binucleated plasma cells, Russell bodies, fibrin, and neutrophils were all significantly increased in the high inflammatory samples (Figure 6A). These findings were consistent with the observed fractions of cells inferred from the gene expression data (Figure 4C). Similarly, the ESR and the levels of CRP, RF, and anti-CCP were increased in patients with high inflammatory scores (Figure 6B).

We also compared the fine specificity of ACPAs between the synovial subtypes (see Supplementary Figure 3, http://onlinelibrary.wiley.com/doi/10.1002/art.40428/abstract), and identified significantly increased levels of autoantibodies to citrullinated filaggrin, vimentin, and fibrinogen a peptides in patients in the high inflammatory subgroup (Figure 6C). Unexpectedly, these subtype differences were not observed when compared according to clinical assessments of pain, tender joint counts, swollen joint counts, or disease duration (Figure 6B). Patients in the low inflammatory subtype had high pain severity scores (median pain score of 6 on a scale of 1–10) but little inflammation in the tissue (according to gene expression data, histologic assessments, or levels of systemic inflammation markers in the blood [i.e., ESR and CRP levels).

We therefore hypothesized that pain might be driven by distinct mechanisms in the various synovial subtypes. To explore this hypothesis, we compared the Spearman’s rank correlation coefficients between the levels of acute-phase reactants and pain scores across all synovial samples grouped together and when parsed according to patient synovial histologic subtype. We found a weak, nonsignificant correlation between pain scores and levels of acute-phase reactants when all patients were analyzed together (Figure 6D), and this became more clear when we divided the patients according to synovial histologic subtype. Patients with high inflammatory synovium, but not those with the other synovial subtypes, showed a stronger, significant correlation between pain scores and CRP levels. This suggests that pain is associated with inflammation in patients with the high inflammatory subtype, and that pain may be driven by distinct mechanisms in the other patients.