Introduction

Autoimmune rheumatic diseases (AIRDs) are systemic diseases characterised by orchestration of disrupted self-tolerance of immune systems, including rheumatoid arthritis (RA) and systemic lupus erythematosus (SLE).1 2 Given shared immunological backgrounds across AIRDs, glucocorticoid and immunosuppressive drugs have been conventionally used for treatment. Recently, new treatments such as molecular-targeted therapy with biological or targeted synthetic disease-modifying antirheumatic drugs (b/tsDMARDs) or JAK inhibitors are successfully indicated for multiple AIRDs.3 However, we face a challenge in characterising clinical immune features and selecting best treatment among multiple targeted therapies for a particular AIRD case towards personalised medicine.4–6 Namely, there exist fundamental but not-yet-solved questions of AIRD aetiology; how should we categorise clinical subtypes of AIRDs, define their heterogeneity within and across AIRDs, and plan effective treatment strategy targeting clinical features of AIRD subtypes? Human genetics of AIRDs provided insights on prominent roles of Polygenic Risk Score (PRS) based on the genome-wide association studies (GWAS).7 8 This achieved unsupervised clustering of AIRDs from genetic backgrounds9 and enabled prediction of a lifetime AIRD risk of individuals.8 However, genetics-based prediction is currently neither modelling underlying clinical heterogeneity within AIRDs, nor reflecting progressing immune profiles in a case.

Immunophenotypes, immune-related cell types in peripheral blood mononuclear cells (PBMCs) quantified by flow-cytometry analysis, can dynamically and quantitatively describe immune profiles of individuals. Flow-cytometry is useful to distinguish: (1) differentiation stage such as naïve T cells or memory T cells, (2) lineage or functional difference such as Th1, Th2 or Th17, (3) activation status or cellular signalling pathways of lymphocytes, which often reflect pathological processes in organs/tissue. Further, immunophenotypes are assayed in a less invasive manner without conflicting the ethical issues of individual identification like genome sequencing. They have advantages in deciphering associations with clinical features of AIRDs such as disease activity and responsiveness to bDMARDs,10–13 being considered as a promising tool towards personalised medicine of AIRD. Despite its utility, however, existing research of immunophenotypes are limited in their sample sizes, coverage of AIRD spectrum and cell types. Further, lack of the standardisation protocol in immunophenotyping assays have made pan-disease or cross-cohort analysis to be difficult.

Here, we introduce our large-scale and cohort-based immune phenotypes dataset that includes over 1000 individuals of a variety of AIRDs and healthy donors attached to deep and longitudinal clinical information. Our aim is to elucidate complex immunonetwork of AIRDs via the multimodal resolutions of immune cells, diseases and individual patients. As principal strategy, immunophenotypes and human genetics are two major components representing the hypothesis-free approaches to elucidate human immunology (online supplemental figure 1). By using a subset of the dataset including 11 AIRDs and healthy controls with immunophenotypes of the 46 cell types, we revealed immune cell type networks associated with AIRDs in a disease-specific manner, as well as shared and distinct immunological features across AIRDs. Unsupervised learning of immunophenotypes successfully dissected the AIRD cases into several clusters, where heterogeneity within a single AIRD, as well as shared architecture across multiple AIRDs, are revealed. By adopting RA as a flagship disease, we found that patient subtypes defined by immunophenotypes could be characterised by different clinical features. Further, immunophenotypes and clinical predisposition of the patients with RA were linked to inborn human genetics.

MethodsStudy participants

We used the data collected through the FLOW study managed by the First Department of Internal Medicine, University of Occupational and Environmental Health, Japan (UOEH). The FLOW study is a prospective cohort recruiting the patients diagnosed with a variety of AIRDs.10–13 In this study, 1088 cases affected with any of the 11 AIRDs registered from May 2011 to December 2018, including RA, SLE, systemic sclerosis (SSc), ANCA-related vasculitis (AAV), idiopathic inflammatory myopathy (IIM), psoriasis, IgG4-related disease (IgG4RD), mixed connective tissue disease (MCTD), ankylosing spondylitis (AS), Sjogren’s syndrome (SjS), and giant cell aortitis (GCA), and 64 controls in the UOEH hospital and the affiliated hospitals were enrolled (table 1). All the patients were diagnosed by professional physicians according to the clinical criteria.14–25 A part of the patients with RA was also registered at the FIRST registry, a prospective cohort longitudinally observing the patients with RA who were newly prescribed with b/tsDMARDs.26 27 The patients affected with the multiple target AIRDs were excluded, while secondary SjS were included considering its relatively high complication rates with the primary AIRDs.

Table 1

Characteristics of the subjects with immunophenotypes

Immunophenotyping

Immunophenotypes of the 46 immune cell types of the participants were longitudinally generated by analysing flow cytometry of PBMCs according to the protocol described elsewhere.10 11 Briefly, after separation and antibody staining (online supplemental table 1), PBMCs were analysed by multicolor flow cytometry (FACSVerse; BD Bioscience). The phenotypes of immune cell subsets were defined based on the Human Immunology Protocol of comprehensive 8-colour flow cytometric analysis proposed by National Institutes of Health/Federation of Clinical Immunology Sciences,28 with minor modifications for detecting Tfh cells. Details of the gating strategy are in online supplemental figure 2. After quality control (QC), immunophenotypes of the 947 patients were obtained.

To normalise potential batch effects among the samples and assays, we developed a pipeline to apply automated counting of the immune cells using R Bioconductor packages of openCyto (version 2.2.0). The list of the 46 cell types is shown in online supplemental table 2. Counts of each cell and sample were log transformed after normalised by the total viable PBMC counts (ie, sum of the CD4+ T cells, CD3-CD19+ B cells, CD3−CD14−CD19−CD20−CD56+ NK cells, CD3−CD14−CD19−CD20−CD56−HLA-DR+ dendric cells (DCs) and CD14+ monocytes). Details of the immunophenotyping including sample collection and immune cell sorting are described in online supplementary methods.

Statistical analysis of the immunophenotypes

Statistical analysis of the immunophenotypes was conducted by using R statistical software (V.4.0.2), as described elsewhere.29 This included Pearson’s correlation analysis, hierarchical clustering implemented in the stats package, principal component analysis (PCA) and logistic regression. Hierarchical clustering of the immune cells, diseases and individuals were conducted based on Pearson’s correlations, Euclidian distances and regression coefficients of the two-dimensional matrix of the normalised immune cell counts respectively, of which axes consisted of the samples and immune cell type categorisations. The matrix was regressed by potential confounding factors and immunophenotyping batch effects before hierarchical clustering, including age, sex and the top 10 PCs. The numbers of the defined clusters were sequentially increased until satisfying the following conditions: (1) clusters could distinguish the major immune cell type classification of T cell, B cell and innate immune cells (for immune cell types), (2) each cluster included at least multiple AIRDs (for diseases), (3) multiple clusters included less than 20% of the patients (for patients). Stability of the clustering results were confirmed by repeating de novo clustering by subsampling the randomly selected 95% of the patients from each AIRD (97.4% concordance of the clustering results by×1000 iterations (online supplemental figure 3). The network graphs were visualised by using the igraph package implemented in R statistical software (V.1.2.6). Statistical significance of the observed associations was adjusted by using Bonferroni’s correction or Benjamini-Hochberg false-discovery rate (BH-FDR).

Clinical and genetic associations with immunophenotype in RA

For the patients with RA registered at both the FLOW study and FIRST registry (n=214), clinical phenotypes including biomarkers and Disease Activity Scores (DAS) were obtained at the time of starting and after 24 weeks of the b/tsDMARDs treatment (see details in online supplemental table 3). Genome-wide genotypes of the patients with RA were obtained by using Infinium Asian Screening Array (Illumina). Details of the genotype calling, genotype QC filters, whole-genome genotype imputation and post imputation QC (minor allele frequency >0.5% and Rsq>0.7) are described elsewhere.30 We estimated PRS of the subjects based on (1) the transethnic RA case–control GWAS (outside of the HLA region at 6p21)8 and (2) the GWAS on RA-interstitial lung disease (ILD),31 using PRS-CSx32 by adjusting potential population stratification using the PCA components.

ResultsCohort-wide Immunophenotyping of the patients with AIRD

Our study enrolled the 1088 patients affected with any of the 11 AIRDs registered at the FLOW study, a prospective cohort with longitudinal follow-up of the clinical information, including RA, SLE, SSc, AAV, IIM, psoriasis, IgG4RD, MCTD, AS, SjS and GCA, as well as the 64 controls (table 1).10 11 We analysed immunophenotypes by multicolor flow cytometry of PBMCs obtained from the participants according to the Human Immunology Protocol (online supplemental table 1 and online supplemental figure 1).28 To normalise potential batch effects, we developed and applied an automated gating and counting pipeline of the immune cells. This advantageously enabled standardisation of immunophenotyping assays across multiple diseases with a cohort-wide scalability. After QC to remove the samples with low quality or technical difficulty in flow cytometry sorting, immunophenotypes were obtained from the 947 subjects as normalised counts of the 46 immune cell types (RA, n=285; SLE, n=170; SSc, n=131; AAV, n=94; IIM, n=77; psoriasis, n=35; IgG4RD, n=35; MCTD, n=22; AS, n=13; SjS, n=20; GCA, n=11; controls, n=54; table 1 and detailed definitions of the immune cell types in online supplemental table 2.

Immune cell type network identified by the cohort-wide immunophenotyping

We first conducted unsupervised hierarchical clustering of the 46 immune cell type abundance defined by the cohort-wide immunophenotyping. The FLOW study recruited the patients with AIRD at initial disease diagnosis or before intensive treatment after disease flare. Therefore, immunophenotypes of the patients with AIRD would reflect primary active states of the disease. We obtained the six major immune cell type clusters. They consisted of the four clusters of T cell subsets (clusters 1, 2, 4, 5), B cell subsets (cluster 3) and innate immune cell subsets (natural killer cells (NKs), DCs and monocytes; cluster 6; figure 1A). The T cell clusters mostly corresponded to those of (1) regulatory T cells (Treg) and follicular helper T (Tfh) cells (cluster 1), (2) activated T cells (cluster 2), (3) CD8+T cells (cluster 4) and (5) CD4+T cells and other T cells (cluster 5). These clusters were concordant to the prior knowledge on peripheral immune cell classifications and their biological kinships,13 supporting robust validity in methodological strategy and data quality of the immunophenotyping. Plasmablasts belonged to the activated T cell cluster differentially from other B cell subsets. This may reflect differentiation induction of B cells into plasmablasts accelerated by activated T cells.

Figure 1

Immunophenotype network across immune cell types. (A) Pairwise correlation of the normalised counts of the 46 immunophenotypes across the patients with AIRD. Cell types were ordered according to the hierarchical clustering analysis. Positive and negative correlations are indicated in red and blue, respectively. Significant correlations satisfied Bonferroni’s multiple testing for the number of the cell type pairs are highlighted (α=0.05). (B) Correlation network of the immunophenotypes. Positive and negative correlations are indicated in red and blue, respectively. Correlations with Pearson’s correlation coefficient |r|>0.4 are indicated. AIRD, autoimmune rheumatic disease; DC, dendritic cell.

Network analysis based on pairwise correlation matrix across immune cell type abundance visualised cell type connections across the patients with AIRD. Tight connection within the major immune cell types of T cells, B cells, DCs, NKs and monocytes were obtained (figure 1B). Again, the network analysis characterised close connection of plasmablasts with activated T cells (eg, activated CD4+, CD8+ and CCR6+CD8+ T cells) rather than with other B cells. We observed negative correlation between CD4+T cells and monocyte lineages, reflecting dominant balance in abundance between these two major cell fractions in peripheral immune cells.

Cross-autoimmune disease network in sight of immunophenotypes

We then conducted unsupervised hierarchical clustering of the 11 AIRDs based on pairwise Euclidian distances across AIRDs calculated from principal components (PCs) of the 46 immune cell types. We obtained the four major AIRD clusters, consisting of SLE and MCTD (cluster 1), AAV, IIM, GCA (cluster 2), controls, SjS, SSc, IgG4RD (cluster 3) and AS, psoriasis, RA (cluster 4; figure 2A). Network analysis visualised connections across the AIRDs. In particular, we observed tight positive correlations (1) between SLE and MCTD and (2) between AAV and IIM, which represented relatively similar immunophenotype features within these AIRD pairs (figure 2B).

Figure 2

Immunophenotype network across AIRDs. (A) Pairwise Euclidian distance across the AIRDs. AIRDs were ordered according to the hierarchical clustering analysis. Euclidian distances are coloured as in the legend. (B) Euclidian distance network of the AIRDs. Network weights were defined as inverses of the Euclidian distances. Network weights more than one-third of the maximum weight are indicated. AAV, ANCA related vasculitis; AIRD, across autoimmune rheumatic disease; AS, ankylosing spondylitis; GCA, giant cell aortitis; IIM, idiopathic inflammatory myopathy; MCTD, mixed connective tissue disease; RA, rheumatoid arthritis; SjS, Sjogren’s syndrome; SLE, systemic lupus erythematosus; SSc, systemic sclerosis.

Immunophenotype network disentangled immune cell type-specificity of AIRDs

To disentangle underlying immunophenotype features defining the AIRD clusters, we constructed immune-phenotype network across the 11 AIRDs and the 46 immune cell types (figure 3). This network visualised immune cell type-specificity of AIRDs, providing immunological explanations of the AIRD clusters obtained in figure 2.

Figure 3

Immunophenotype network across AIRDs and immune cell types. AIRD-immune cell type associations based on backward-forward stepwise logistic regression adjusted by sex were visualised as a network plot. Associations that satisfied Bonferroni’s correction (α=0.05), BH-FDR (<5%, 10%), and nominal significance (p<0.05), were shown as edges with the corresponding colours as indicated in the legend. Dotted and solid lines represent decreasing and increasing directional dosages of the immune cells on AIRDs, respectively. AAV, ANCA related vasculitis; AIRD, across autoimmune rheumatic disease; AS, ankylosing spondylitis; BH-FDR, Benjamini-Hochberg false-discovery rate; SLE, systemic lupus erythematosus; SSc, systemic sclerosis.

MCTD is characterised as Raynaud’s phenomenon, anti-U1-RNP autoantibody and multiple symptoms of AIRDs,33 and the obtained results indicate tight connections in the network with SLE. Enhancement of activated T cell subsets characterised these two diseases. Contrarily, SSc belonged to the same cluster as IgG4RD and SjS. IgG4RD is a novel disease entity characterised by serum IgG4 elevation and tissue infiltration by IgG4-positive plasma cells, followed by symptoms related to multiple AIRDs including SjS and SSc.34 35 Increase of effector memory CD4+T cells, Tfh cells, and plasmablasts and decrease of naïve CD4+T cells distinguished this cluster. Clinical manifestations of AAV and IIM are diverse,36 so their relatively shared immunophenotype was out of expectation. Increase of B cells and monocytes were commonly observed. RA, psoriasis, AS were clustered together, and characterised with autoimmune joint arthritis and shared clinical features such as indication of anti-TNF therapy.37 38 These diseases demonstrated as decrease of activated T cells, Treg, Tfh cells, plasmablasts and monocytes. Abnormality of neutrophil extracellular traps has been reported to play a role in the pathogenesis of SLE and AAV.36 39 40 We note that this study did not investigate neutrophils, which are not included in PBMCs.

In summary, immunophenotype-based clustering reflected shared clinical features within the defined AIRD clusters, as well as their heterogeneity across the clusters. Underlying immune cell type-specificity prescribes the clusters’ characteristics.

Patient-level clustering of AIRDs based on immunophenotypes

Motivated by successful de novo classifications of AIRDs based on immune cell type specificity, we then conducted unsupervised clustering of individual patients from a cohort-wide immunophenotypes without prior clinical information other than age and sex. Hierarchical clustering based on PCs obtained the major six patient clusters (figure 4A). We assessed patient characteristics separately for each AIRD and found conspicuous heterogeneity in the patient distributions across the clusters for most of the AIRDs (p<0.05/12=0.0042 for MCTD, SLE, AAV, IIM, control, RA; p<0.05 for SSc, AS and psoriasis). Hierarchical clustering of AIRDs based on the patient-level distributions again demonstrated similar AIRD classifications to those obtained in the disease-level clustering in figure 2A (eg, tight connection between AAV and IIM and among AS, psoriasis and RA; figure 4B). Thus, patient-driven clusters of AIRDs could reflect underlying immune cell type specificity of AIRDs, as well as their shared and distinct biological backgrounds, which were identified through our analysis mentioned in the previous sections.

Figure 4

Patient-level AIRD clustering. (A) Unsupervised clustering of individual patients with AIRD according to the immunophenotypes. (B) Clustering of AIRD based on patient-level clustering results. AIRDs with significantly heterogeneous patient distributions across the clusters were highlighted. †Bonferroni’s correction (α=0.05) and *p<0.05. (C) Pie charts showing proportions of the AIRDs within each cluster, and those of the clusters within RA and SLE. (D) A co-plot of the fractions of the patients with RA and SLE within each cluster. AAV, ANCA-related vasculitis; AIRD, across autoimmune rheumatic disease; AS, ankylosing spondylitis; GCA, giant cell aortitis; IIM, idiopathic inflammatory myopathy; MCTD, mixed connective tissue disease; RA, rheumatoid arthritis; SjS, Sjogren’s syndrome; SLE, systemic lupus erythematosus; SSc, systemic sclerosis.

We then assessed fractions of the AIRDs to each patient cluster. We particularly focused on RA and SLE (n=285 and 170, respectively), the two representative AIRDs which consist of the largest sample sizes in our immunophenotype registry with significantly heterogeneous distributions of the patients across the clusters (figure 4C). We found that the clusters 1–3 showed higher fraction of the patients with SLE than in the patients with RA. In contract, the clusters 4–6 showed higher fractions in RA than in SLE (figure 4D). This briefly projects general differentiation of the patients with AIRD to those with RA- or SLE-like immunophenotypes. Namely, our data-driven analysis empirically confirmed the well-recognised clinical observations that RA and SLE are very different diseases.1 2 Of note, relatively small proportions of the patients with RA or SLE (23.2% or 18.8%) belonged to the non-RA or non-SLE-like clusters. Discussions on patient heterogeneity and treatment strategies in RA or SLE have been long-standing.41 This result provides a piece of empirical evidence of disease subtypes in a single AIRD from the viewpoint of the patient-oriented immunophenotypes.

Immune cell type and clinical manifestations of the RA patient clusters

To examine how individual patients’ immunophenotypes define clinical characteristics of the heterogeneous AIRD subsets, we conducted association studies between clinical information and patient clusters. We focused on the patients with RA, the AIRD with the largest sample size in our datasets. Further, parallel registration of the patients with RA in another clinical cohort of the FIRST registry enabled us to investigate detailed longitudinal follow-up including disease activity and treatment responses (online supplemental table 3).26 27

Each cluster demonstrated specific clinical features (figure 5, online supplemental tables 3 and 4). We did not observe apparent differences of disease durations and female-male ratios across the clusters (FDR>0.05), or treatment choices of the four major biologics categories (TNFα inhibitors, JAK inhibitors, CTLA4-Ig and IL6R inhibitors; p>0.23). As for the patients with RA in the cluster 1, the largest AIRD patient cluster with SLE-like immune-phenotypes, were associated with decreased activated Tregs, Tfh cells, plasmablasts and Th17 cells, and increased Th1 cells and CD4+T cells. The cluster 1 showed average disease activity at enrolment before treatment. However, this cluster was associated with poor treatment response in Evaluator’s global assessment (EGA) after 24 weeks of treatment, showing treatment resistance. Majority of the patients with RA belonged to the cluster 4, the cluster with RA-like immune-phenotypes. The patients with RA in this cluster were characterised as high methotrexate dosage at baseline but showed improvement in EGA and morning stiffness after 24 weeks, presenting its initially severe disease activity but better treatment response. Associated immunophenotypes of the cluster 4 were similar to those of the cluster 1, but negative associations with the T regs were specific to the cluster 1. This may explain treatment resistant feature of the cluster 1, as well as its SLE-like characteristics. We note that no differences in anticitrullinated protein antibody (ACPA) and rheumatoid factor (RF) positivity across the clusters (on average 74.3% and 90.0%, respectively). ACPA positivity is known as a clinical biomarker to predict joint destruction progression and organ failures.42 Immunophenotypes could reflect clinical heterogeneity in disease activity as well as treatment response beyond these classical biomarkers. As for the patients with SLE, We observed significant differences of the renal diseases across the SLE patient clusters (FDR<0.005). Lower renal diseases were observed in the clusters 4 and 5 where RA-like immunophenotypes were dominant (16.0%), when compared with the other clusters (46.2%) (online supplemental table 5).

Figure 5

Immune cell type and clinical information associations with the RA patient clusters. Immune cell types and clinical information with significant associations with the RA patient clusters are indicated. Outer and inner circles represent the AIRD and RA patient sample sizes belonging to the clusters, respectively. Edges represent significance levels and directional associations as in the legend. The Δ values show improvement in the disease activity values after treatment (ie, larger positive Δ values corresponding to better treatment responses). The list and values of the clinical information is available in online supplemental tables 3 and 4. AIRD, across autoimmune rheumatic disease; CRP, C reactive protein; DAS-28, Disease Activity Score-28; RA, rheumatoid arthritis.

Human genetics involvements in clinical and immunological manifestations of RA

The aetiopathogenesis of the AIRDs is in part human genetics,29 43 which is considered as another promising tool of data-driven approach to decipher disease backgrounds. However, how genetics affect clinical and immunological manifestations within the patients with AIRD after disease onset has been elusive. To investigate human genetics involvements in clinical and immunological manifestations of RA, we conducted an association study of PRS, a genome-wide aggregation of disease risk variants of the individuals to predict phenotype risk,44 45 with clinical information and immunophenotypes of the patients with RA registered in the FIRST registry26 27 (n=123; figure 6).

Figure 6

Associations of PRS with clinical and immunological information. Associations of the PRS based on RA case-control GWAS8 and within RA case GWAS for RA-ILD31 with clinical information and immunophenotypes of the cell types are indicated. Red and blue colours represent positive and negative associations as in the legend, respectively. FDR, false-discovery rate; GWAS, genome-wide association studies; ILD, interstitial lung disease; RA, rheumatoid arthritis; PRS, Polygenic Risk Score.

We first assessed PRS based on the latest cross-population RA case-control GWAS,8 a promising surrogate to predict RA onset susceptibility. RA case–control PRS was associated with younger age of RA onset, autoantibody positivity (ACPA and RF), higher disease activity (DAS-C reactive protein, DAS-ESR, Simple Disease Activity Index and DAI), and poorer treatment responses after 24 weeks (BH-FDR<0.10). Negative associations with decrease of naïve Treg were observed (p<0.05). Thus, PRS based on RA disease onset risk could also be useful to predict disease activity and prognosis.

We then assessed PRS based on the stratified GWAS within the RA cases. We adopted the Japanese GWAS of RA-ILD,31 one of the clinically important complications with high morbidity and mortality.46 While RA-ILD PRS did not significantly associate with clinical manifestations, we observed positive associations with increased DCs (myeloid DCs, DCs and plasmacytoid DCs; BH-FDR<0.10). While detailed aetiology or RA-ILD has been fully investigated, involvement of DCs in ILD in the rheumatic disease model have been reported.47 Our results provided additional supportive evidence in the human RA patients.

Discussions

In this study, we evaluated the large-scale and cohort-based immunophenotype dataset over 1000 individuals and 11 AIRDs with deep and longitudinal clinical information. De novo clustering of immune cell types and AIRDs successfully constructed classifications consistent with prior knowledge on immunological features of AIRDs. We obtained immune cell type spectra specific to each AIRD or shared across different AIRDs. Further, this provided novel insights into disease aetiology and their underlying networks as well. As an insightful observation, we observed close network between MCTD and SLE. In contrast, SSc was clustered together with SjS and IgG4RD, AIRDs characterised with tissue fibrosis and autoantibody production. This may motivate potential utility of the treatment targeting B cell lineages for such AIRDs. Nevertheless, our findings are based on peripheral immune cells, and further investigation using the lesion tissue is warranted. Detailed immunocharacterisation specifically focusing on each AIRD and longitudinal clinical information would also be a next step. While our cohort focuses on primary status of the AIRDs, we admit potential impact of treatment choices and prescription behaviour on immunophenotypes of the patient. Sex-dependent differences in immunological and clinical features are also of interests. While we showed stability, further validation of the hierarchical clustering results would be warranted.

Individual patient-level clustering simplified the patients with AIRD into several clusters, where RA-like and SLE-like ones were exclusively dominant. The immunological abnormalities found in the patients with RA and SLE reflect their clinical differences (eg, incidence of SLE-like symptoms in patients with RA as consequence of anti-TNF therapy). Since we focused on primary disease status of AIRD and did not handle secondary SjS as an independent disease, elucidation of immune profile shift driven by secondary SjS would be also warranted. In general, patients are likely to manifest more than one autoimmune condition. Such complex conditions may make patient diagnosis based on immunophenotyping challenging.48

In-depth analysis of the patients with RA revealed that such patient clusters differentially defined clinical heterogeneity in disease activity and longitudinal treatment responses. Of interest, the patients with RA with SLE-like immunophenotypes such as decreased Treg were characterised as treatment resistance. We note that the FLOW study mostly focused on recruitment of the early-onset patients with RA with high disease activity before treatment (bio-naïve for >75% of the patients). This mitigated biases from patient heterogeneity in disease status or treatment strategies and empowered the analysis differently from previous studies. Our results had advantageous messages that cohort-wide and cross-disease immune-phenotyping would be required to pinpoint such heterogeneous patient subtypes existing within the single disease. Given that our approach identified the disease subtypes along with immune cell type specificity, this facilitated personalised treatment strategy or novel drug discovery targeting such key immune cell types.49

Inborn human genetics highlighted with PRS showed associations with dynamics in clinical manifestations and immunophenotypes. Previous studies reported values of PRS to predict disease onset in early stages of life,44 45 but how PRS could benefit patient stratification and personalised treatment strategy after onset is not yet widely validated. We showed that PRS from both case–control GWAS and within-case stratified GWAS could contributed to fulfil such missing pieces. Further, associations between GWAS-derived PRS and immunophenotypes suggested candidate key immune cell types implicated in disease aetiology (eg, DCs for RA-ILD in figure 6), proposing a novel biological utility of PRS.

While our study achieved translayer integrations across immunophenotypes, clinical information and genetics, further involvement of multilayer omics resources should be warranted to more finely characterise AIRDs. This may include metabolome,50 proteome,51 microbiome,52 virome53 and somatic alterations.44 One potential limitation of immunophenotyping is being unable to know transcriptomic and epigenetic dynamics within the immune cells. Recent development of single cell technologies could complement this layer, and its application in the field of AIRDs are expected.54 Nevertheless, expansion of the target diseases into wider ranges of AIRDs and immune-related human complex phenotypes is also the future step. Dynamic and longitudinal effects of treatment on immunophenotypes is of clinical interests. We note that the number of the control subjects (n=64) were relatively smaller than those of the AIRDs. Enlargement of the control subjects is desirable to further validate our findings on AIRDs. Finally, while our study provided robust associations of immunophenotype with AIRD clinical features, we need more evidence to reliably predict patients’ future clinical outcomes.

In summary, our study provided a piece of evidence that large-scale immunophenotypes could characterise AIRDs, define underlying within-disease subsets, and profile individual patients, which would be useful for implementing personalised medicine of AIRDs.