Genomic risk scores for juvenile idiopathic arthritis and its subtypes

Training and external validation of the JIA genomic risk score

We used 10×10-fold cross-validation to tune the penalised model and estimate the area under the receiver operating characteristic curve (AUC) in the UK dataset, achieving a maximum AUC=0.671 (95% CI 0.668−0.674) and selecting 26 SNPs in the model (Figure S2). The selected SNPs and their weights in the GWAS clearly showed the strong HLA dependency of JIA (Figure S3). The final model (SNP weights) was held fixed for external validation GRS (http://www.pgscatalog.org/pgs/PGS000114/).

We applied the GRS to the CLARITY and CHOP datasets and evaluated the prediction power of our model in terms of AUC and odds ratios (Table 2). Overall, the GRS showed highly consistent performance across both validation datasets, with AUCs of 0.657 and 0.671 for CHOP and CLARITY, respectively. The corresponding odds ratios per s.d. of GRS were 1.831 (95% CI 1.685-1.991) and 2.008 (95% CI 1.731-2.345) in CHOP and CLARITY, respectively. The associations of the GRS with JIA status were unattenuated when adjusting for the top 10 genetic principal components (PCs) and sex. Furthermore, we compared the GRS model with a model generated using known SNPs highly association with JIA (Table 1 from Hinks et al.¹⁶). Using this model, we achieved AUC of 0.614, 0.642 and 0.648 for UK, CHOP, and CLARITY respectively, showing that the GRS model was a stronger risk predictor.

Recent work has also shown that a metaGRS approach can substantially improve genomic prediction of common diseases²¹, thus we compared the metaGRS method to the lasso-penalised regression above (Supplementary Material). We did not observe that a JIA metaGRS, based on a set of related autoimmune disease GWAS summary statistics, significantly improved upon the JIA GRS above, thus subsequent analyses utilised lasso-regression.

Subtype analysis

We next extended our analysis to consider subtypes of JIA and construct subtype-specific GRSs thereof. The International League of Associations of Rheumatology (ILAR) classification system recognises seven subtypes of JIA, systemic arthritis, oligoarthritis, rheumatoid-factor-positive polyarthritis (RF-Positive), rheumatoid-factor-negative polyarthritis (RF-Negative), enthesitis-related arthritis (ERA), psoriatic arthritis, and undifferentiated arthritis⁵. Subtypes vary substantially in average age at onset, sex distribution, number of joint affected, and clinical features (see Figure 2)²⁶. Of particular interest clinically is the time between onset of symptoms to visiting a paediatric rheumatologist which has been estimated to vary from 11 months (range: 8 to 70 weeks) in the case of enthesitis-related JIA to 1 month (range: 2 to 36 weeks) for systemic JIA^14,27. The heterogeneity of JIA was reflected in our case data; all seven subtypes were present, albeit at frequencies ranging from common (∼41% of UK cases were oligoarthritis JIA) to relatively rare (∼4% of UK cases were undifferentiated JIA) (Table 3).

For each of the seven subtypes, we utilised the UK dataset to train subtype specific GRSs, employing a similar lasso-regression in cross-validation approach as for the JIA GRS above (see Methods). Each subtype GRS was trained on that subtype’s cases and all non-JIA controls, with all other JIA subtypes were excluded from training, and then externally validated in CHOP and CLARITY. We found that there was a high degree of variability in discrimination between subtype GRSs, with some subtypes displaying cross-validated AUCs greater than the JIA GRS and others not being significantly different than random chance (AUC=0.5; Figure 2).

• Download figure
• Open in new tab

Figure 2:

10×10 cross-validated AUC (LOESS-smoothed) AUC achieved by training the seven JIA subtype specific models (top) and median time taken by an individual with JIA to be referred for first time to a paediatric rheumatologist visit (in months; bottom)²⁷.

When compared to a previous estimate of the median time for a child with JIA to be referred to a paediatric rheumatologist visit²⁷, we found that the JIA subtype (enthesitis-related) with the longest time-to-referral (∼11 months) also had the strongest subtype GRS (AUC=0.80 with 149 SNPs). This result may be partially due to circulating HLA-B27 antigen being a component of the enthesitis-related subtype classification criteria^5,28. To test this, we utilised previously reported HLA-B27 tag SNPs^29,30 (Supplementary Material), and found that our enthesitis-related JIA GRS still outperformed the HLA-B27 tag SNPs model in the UK dataset (AUC=0.75). The next strongest subtype GRS (AUC=0.71 with 25 SNPs) was for oligoarthritis JIA, which has a median time-to-referral of ∼3 months. While JIA subtype rates are not well characterised, oligoarthritis JIA is the most common clinical subtype reported³¹.

The weakest subtype GRSs were for the undifferentiated (AUC=0.573 with 119 SNPs) and systemic (AUC=0.546 with 3,430 SNPs) subtypes. This was not unexpected as these subtypes are somewhat different to the other five subtypes. Systemic JIA is considered an autoinflammatory disease and, while there is little genetic overlap with the other JIA subtypes³², it has strong genetic signals from the HLA class II molecule encoded by HLA-DRB1*11, confirming the role of the adaptive immune system³³. Children are diagnosed with the undifferentiated subtype when their symptoms do not fit within other subtypes, or meet the criteria for multiple subtypes.

External validation of the subtype GRSs in CLARITY showed highly consistent AUC estimates with cross-validation performance in the UK (Table 4). The CHOP dataset showed somewhat less consistent external validation than CLARITY for subtype GRSs, but in both CLARITY and CHOP the top performing subtype GRSs were both the enthesitis-related JIA GRS (AUCs of 0.83 and 0.69, respectively) and oligoarthritis JIA GRS (0.75 and 0.77, respectively). In both external validations, the systemic JIA subtype GRS had poor predictive capacity with AUCs of 0.47 and 0.55 in CLARITY and CHOP, respectively.

Ethics Statement

All participants gave informed consent and the study protocols were approved by the relevant institutional or national ethics committees. The Australian CLARITY cohort collection was approved by the Royal Children’s Hospital Human Research Ethics Committee; UK ethical approval was obtained from the North West Multi-centre for Research Ethics Committee (MREC:02/8/104 and MREC:99/8/84), West Midlands Multi-centre Research Ethics Committee (MREC:02/7/106), North West Research Ethics Committee (REC:09/H1008/137), and the NHS Research Ethics Committee (REC:05/Q0508/95); and the US CHOP cohort collection was approved by the institutional review boards of the Texas Scottish Rite Hospital for Children, the Children’s Mercy Hospitals and Clinics, and the Children’s Hospital of Philadelphia.

Phenotype and clinical data

The International League of Associations for Rheumatology (ILAR) classification system⁵ provides generally accepted guidelines for researchers and clinicians to delineate the seven mutually-exclusive categories of JIA based on the predominate clinical and laboratory features. For example, to diagnose the enthesitis-related JIA subtype, an individual has have at least two of the following conditions: presence of HLA-B27 antigen; onset of arthritis in a male over age six; acute anterior uveitis; presence of or history of sacroiliac joint tenderness and/or inflammatory lumbosacral pain; or have a first-degree relative with history of ankylosing spondylitis, enthesitis-related arthritis, Reiter’s syndrome, or acute anterior uveitis.

In this work, the JIA diagnosis for the three cohorts was made according to the ILAR revised criteria and all JIA cases were of age of onset <16 years old. In the UK cohort¹⁶, the JIA cases were obtained from five sources: The British Society for Paediatric and Adolescent Rheumatology National Repository of JIA; a group of UK cases with long-standing JIA described previously³⁶; a cohort collected as part of the Childhood Arthritis Prospective Study¹⁰; a cohort of children recruited for the SPARKS-CHARM (Childhood Arthritis Response to Medication)³⁷; and an ongoing collection of UK cases from the UK JIA Genetics Consortium. The UK controls were from the shared UK 1958 Birth cohort and UK Blood Services Common Controls, established as part of the WTCCC³⁸.

For the CHOP cohort^18,39, the clinical data of JIA cases were collected from the electronic health records (EHR) completed by the pediatric rheumatology specialist within the Division of Rheumatology at CHOP and abstracted into a JIA Registry maintained within the Center for Applied Genomics (CAG) at CHOP. Control subjects were unrelated and disease-free children recruited by the CAG recruitment team within the CHOP Healthcare Network. In addition, control subjects had no history of JIA or other chronic illnesses and were screened as negative for a diagnosis of autoimmune diseases, based on data from CHOP’s electronic health record and by intake questionnaires obtained by the recruiting staff from the CHOP Center for Applied Genomics.

Finally, in the CLARITY cohorts^13,40, all the cases were aged 0–18 years at interview, with diagnosis of JIA by a paediatric rheumatologist. Incident cases were defined as children recruited within six months of diagnosis and prevalent cases were defined as those children diagnosed more than 6 months before recruitment and since 1997. Controls were recruited through the Royal Children’s Hospital Day Surgery Unit. Exclusion criteria for cases and controls were the presence of major congenital abnormalities or illness that would forgo school attendance in the one year prior to recruitment.

Genotype data and quality control

All collected cases are pediatric rheumatologist-diagnosed and controls have no history of any autoimmune disease. Additionally, all participants are white European descent with the exception of ∼20% of the participants in CLARITY, who are of non-European ancestries. Finally, all genotyped included in each dataset were aligned to the GRCh37/hg19 assembly build and passed stringent quality control measures (detailed below).

The initial UK dataset consists of 2,758 cases and 5,187 controls. Controls were obtained from the Wellcome Trust Case Control Consortium (WTCCC) which have been demonstrated to be well-matched to UK JIA cases^16,41. UK cases were collected through five cohort studies recruiting from multiple tertiary centers around the UK, where 1,670 were genotyped on the Illumina HumanOmniExpress array and 1,088 on the Illumina HumanCoreExome array. The CHOP (US) dataset consisted of 1,229 cases from US and Norway, and 5,512 pediatric controls, all recruited from within the CHOP Healthcare Network. All CHOP controls and 1,008 JIA cases were genotyped using Illumina HumanHap550 or Human610-Quad arrays⁴⁰ and the remaining 250 cases were genotyped using the Illumina HumanCore array. Lastly, the Australian CLARITY (ChiLdhood Arthritis Risk factor Identification sTudY)¹³ consisted of 558 cases and 704 controls collected from the Royal Children’s Hospital and Monash Medical Centre in Melbourne. All controls and 406 cases were genotyped using the Illumina HumanCore array, with the remaining 152 cases genotyped on the Illumina HumanHap550^17,18.

Before genotype imputation, we applied stringent quality control (QC) procedures standardised across all datasets. The initial QC performed, using plink1.9^42,43, consisted of removing non-autosomal SNPs, SNPs with MAF<1%, SNPs and individuals with missingness>10%, and SNPs with deviations from Hardy-Weinberg Equilibrium in controls (P<10⁻³). Additionally, using KING v2.1.5⁴⁴, we identified and removed one of two individuals with a second or higher degree relatedness within the datasets. The resulting datasets consisted of 7,505 individuals (2,324 cases and 5,181 controls) with 144,964 SNPs for UK, and 6,741 individuals (1,229 cases and 5,512 controls) with 489,672 SNPs for CHOP. The CLARITY dataset consisted of three batches; after quality control, these consisted of 152 individuals (all cases) with 490,629 SNPs, 838 individuals (247 cases and 591 controls) with 246,814 SNPs, and 118 individuals (5 cases and 113 controls) with 243,662 SNPs, respectively.

For genotype imputation, we used the Michigan Server to impute the datasets using Minimac3 and the 2016 Haplotype Reference Consortium (HRC) as the reference panel⁴⁵. We merged all the subsets of CLARITY into a single dataset. Then, for each dataset (UK, CHOP and CLARITY), we removed multi-allelic and duplicated SNPs, SNPs with imputation r²<0.5 and MAF<0.01, SNPs deviating from the Hardy-Weinberg Equilibrium in controls (P<10⁻³), and those with ambiguous strand (A/T or C/G alleles).

Next, we performed principal component analysis (PCA) using FlashPCA2^46,47 over the filtered data (Figures S8–S11). For each dataset we selected the largest homogeneous subset of individuals based on the first five principal components within each set, where we defined as outliers the set of cases (controls) individuals that clearly differed in the PCA plots from the controls (cases). The UK dataset was kept in its entirety; for CLARITY n=168 individuals were removed and for CHOP n=3,228 individuals were removed (including non-US individuals). Figures S12 and S13 show the principal component analysis for the resulting subsets. Table 1 shows the final datasets; for the analysis we used the n=5,545,761 genotyped and imputed SNPs which were available post-QC on all three datasets.

Before training our model we also estimated the heritability and SNP-heritability of JIA in our datasets. The sibling recurrence risk of JIA has previously been estimated to be 11.6-fold¹⁹, and for first cousins 5.8-fold; assuming the JIA prevalence to be 1 in 1000, we estimated JIA heritability to be approximately 54%. For SNP-heritability, we used GCTA version 1.91.7^27,48 and adjusted for the first ten principal components together with an assumed JIA population prevalence of 1/1000. The estimated SNP-heritability was 0.25 (std err 0.02) for UK, 0.37 (std err 0.13) for CLARITY, and 0.51 (std err 0.07) for CHOP.

Development and validation of the genomic risk score

Because it was the most homogenous dataset with the largest case sample size (2,324 cases, 5,181 controls), we utilised the UK dataset to train our models. To account for potential confounding by the case/control genotyping batch in the UK dataset, we used logistic regression of case/control status on sex and the first 10 genetic PCs. The PCs were computed over a subset of the SNPs of UK, excluding the HLA region as well as known or putative JIA risk loci¹⁶ (defined here as SNPs with P<10⁻⁵ and all SNPs within 1 Mb of the former). The residuals from the regression were then used as the phenotype for constructing the GRS.

To create the GRS we used SparSNP²⁵, which is an efficient implementation of a lasso-penalised linear model which has been previously shown to outperform other methods when there are known to be strong effects within regions of complex LD, such as the MHC⁴⁹. This model is induced by minimising the L1-penalised squared loss function (standard lasso) over N individuals and m SNPs, where x_i is the m-vector of genotypes for the ith sample in allele-dosage coding {0,1,2}, y_i is the phenotype, is the m-vector of weights, is the intercept, and λ is the L1 penalty. All SNPs in the selected training dataset are considered in the model but typically only a small number of SNPs receive a non-zero weight. The number of SNPs with non-zero weight varies depending on the value of the L1 penalty (higher penalties lead to fewer SNPs) which was tuned via 10-fold cross-validation. The optimal number of SNPs selected (non-zero weighted) in the chosen model is decided based on the model with the highest average AUC (or R² if the phenotype is continuous) across the replications. We also explored adding an L2 penalty to the model (elastic-net); however, no significant improvement was found.

Once a model was chosen, we computed genomic risk scores (GRS) for each of the test datasets (CHOP and CLARITY). Assuming that the number of SNPs is m, then the GRS g_i for a new individual with genotypes x_i is where and are the SNP weights obtained from the model. Subsequently, in order to validate the GRS, we evaluated it using logistic regression in the CHOP and CLARITY datasets, adjusting for sex and ten first genetic principal components of each test set.