|Home | About | Journals | Submit | Contact Us | Français|
To discover novel RA risk loci, we systematically examined 370 SNPs from 179 independent loci with p<0.001 in a published meta-analysis of RA GWAS of 3,393 cases and 12,462 controls1. We used GRAIL2, a computational method that applies statistical text mining to PubMed abstracts, to score these 179 loci for functional relationships to genes in 16 established RA disease loci1,3-11. We identified 22 loci with a significant degree of functional connectivity. We genotyped 22 representative SNPs in an independent set of 7,957 cases and 11,958 matched controls. Three validate convincingly: CD2/CD58 (rs11586238, p=1×10−6 replication, p=1×10−9 overall), and CD28 (rs1980422, p=5×10−6 replication, p=1×10−9 overall), PRDM1 (rs548234, p=1×10−5 replication, p=2×10−8 overall). An additional four replicate (p<0.0023): TAGAP (rs394581, p=0.0002 replication, p=4×10−7 overall), PTPRC (rs10919563, p=0.0003 replication, p=7×10−7 overall), TRAF6/RAG1 (rs540386, p=0.0008 replication, p=4×10−6 overall), and FCGR2A (rs12746613, p=0.0022 replication, p=2×10−5 overall). Many of these loci are also associated to other immunologic diseases.
Rheumatoid arthritis is a chronic autoimmune disease characterized by an inflammatory polyarthritis12. Genetic studies have now identified multiple risk alleles for autoantibody positive RA within the MHC region, a PTPN22 missense allele, and risk alleles in 14 other loci (see Table 1)1,3-11. Most RA risk loci contain multiple genes, and currently the causal genes are unknown. However, most contain at least one plausible biological candidate gene involved in immune regulation, and these genes suggest an important set of processes involved in RA pathogenesis. For example, risk alleles highlight genes involved in T-cell activation by antigen presenting cells (class II MHC region, PTPN22, STAT4, and CTLA4), the NF-κB signaling pathway (CD40, TRAF1, TNFSF14, and TNFAIP3, and the recent report of REL13), citrullination (PADI4), natural killer cells (CD244), and chemotaxis (CCL21).
Based on these observations, we hypothesized that as yet undiscovered autoantibody positive RA risk loci might also contain genes with functions similar to those of genes in known risk loci. Therefore, known RA risk loci can be used to prioritize SNPs for replication from GWAS, especially SNPs with modest statistical support, in independent samples (Figure 1).
To objectively quantify the degree of functional similarity between genes within candidate loci (identified from GWAS) and genes within validated RA risk loci, we used a published functional genomics method, GRAIL (Gene Relationships Across Implicated Loci)2. GRAIL quantifies functional similarity between genes by applying established statistical text mining methods14 to text from a reference database of 250,000 published scientific abstracts about human and model organism genes. For each locus, GRAIL identifies the gene with the greatest number of observed relationships. GRAIL estimates the statistical significance of the number of observed relationships with a null model where relationships between genes near SNPs occur by random chance. This significance score, ptext, represents the output GRAIL score. GRAIL is already able to effectively identify functional inter-connectivity between genes within the previously known RA loci (Figure 2); it might also be able to establish connections between these 16 loci and as yet undiscovered RA risk loci.
Since GRAIL might demonstrate variable performance across different phenotypes, we wanted to carefully quantify its predictive ability in RA before using it to prioritize SNPs for replication. To estimate GRAIL’s ability to distinguish true RA loci from spurious associations, we examined 12 RA risk loci discovered since 2006 (see Table 1, Supplementary Table 1 online). The current GRAIL implementation is based on PubMed abstracts published prior to December 2006. As these 12 risk loci were subsequently discovered, they constitute a representative set to evaluate GRAIL’s performance. In a leave-one-out analysis, we used GRAIL to score each of these loci for functional relationships to the other 15 validated RA risk loci. A total of 10 of the 12 loci obtain GRAIL scores of ptext<0.01. This analysis suggests that at that ptext threshold, GRAIL has an ~83% true positive rate (or sensitivity). To assess the false positive rate of this same ptext threshold, we modeled spurious loci by sampling 10,000 random SNPs from the Affymetrix 500K array; we scored these SNPs against all 16 RA loci. Of the random SNPs, 5.4% scored ptext<0.01; this corresponds to a specificity of ~95%. Assessment of true and false positive rates at different cutoffs revealed an area under the curve (or C statistic) of 0.97 (see Supplementary Figure 1). We note that if a large number of candidate SNPs are screened in a study, this might still result in a large number of false positives.
Next, we attempted to identify novel RA risk loci from a set of SNPs with modest evidence of association from our recent GWAS meta-analysis of 3,392 cases and 12,462 controls1. In our original study, we genotyped SNPs with p<10−4 in the meta-analysis, and found that 6 out of 31 replicated in our independent samples. However, many RA risk alleles have modest effects (e.g., OR <1.2) and will be missed at that significance threshold. We therefore expected that some SNPs at p<0.001 may be risk alleles. After excluding SNPs that were known validated RA risk loci, we identified a total of 370 SNPs from 179 independent regions that obtained p<0.001 (see Methods and Supplementary Note online). The total number of SNPs observed at this threshold is consistent with the approximate number of SNPs expected by chance, suggesting that the majority represent spurious associations and should not be reproducible in an independent case-control study.
For each of the 179 candidate loci we selected the single SNP with the strongest evidence from the GWAS meta-analysis, and then scored it against the 16 validated RA risk loci with GRAIL. If all 179 SNPs were spurious, then ~10 should score ptext<0.01 based on the estimated 5.4% false positive rate. However 22 of the 179 (12.3%) scored ptext<0.01 (see Figure 3A, Supplementary Table 2 online). This represented a significant enrichment compared to random sets of 179 SNPs (p=3.3 × 10−4 by simulation). We therefore expected that of this select subset of 22 SNPs, as many as half might represent true RA risk alleles.
In order to identify which of these 22 SNPs represent true RA risk loci, we genotyped them in an independent validation study of 7,957 cases and 11,958 controls from 11 collections from Europe and North America (see Supplementary Table 3). All cases met 1987 American College of Rheumatology classification criteria15 or were diagnosed by a board-certified rheumatologist, and were seropositive for disease-specific autoantibodies (anti-cyclic citrullinated peptide [CCP] antibody or rheumatoid factor [RF]). All individuals were self-described white and of European ancestry. We assessed association with a Cochran-Mantel-Haenszel (CMH) stratified association statistic16. For each SNP we calculated a z-score, where a z>0 indicates the same allele confers risk in both the replication and the meta-analysis samples. To interpret statistical significance, we used a Bonferroni-corrected one-tailed p-value of 0.0023 (=0.05/22, z>2.83). Additionally, we calculated the overall association p-value across all samples (GWAS meta-analysis plus replication).
Strikingly, of the 22 SNPs examined, 19 (86%) obtained z>0 (see Figure 3B). If these SNPs represented spurious associations then only about half should have z>0; the probability of such a positive skew by chance alone is pskew=0.0005 − suggesting the likelihood of a large number of true RA risk loci within this set of 22 SNPs.
Of the 22 SNPs selected by GRAIL, 13 obtained nominal levels of association at p<0.05 (corresponding to z>1.65); whereas no more than 2 might be expected by chance alone. More compellingly, 7 SNPs achieved a Bonferroni-corrected level of significance in replication (p<0.0023, z>2.83).
When we aggregated both GWAS meta-analysis and replication genotype data (see Table 2, Supplementary Table 4 online), we observed the strongest associations at rs11586238 on 1p13.1 near the CD2 and CD58 genes (1.4×10−6 replication, 1.0×10−9 overall), at rs1980422 on 2q33.2 near CD28 (4.7×10−6 replication, 1.3×10−9 overall), and at rs548234 on 6q21 near PRDM1 (1.2×10−5 replication, 2.1×10−8 overall). Based on conservative estimates of genome-wide significance these represent confirmed RA risk alleles.
Four additional loci replicated; however, aggregate analysis of GWAS meta-analysis and replication genotype data did not exceed a conservative estimate of significance. We observed evidence of association at rs394581 on 6q25.3 near TAGAP (1.5×10−4 replication, 3.8×10−7 overall), rs10919563 on 1q31.3 within a PTPRC intron (2.6×10−4 replication, 6.7×10−7 overall), rs540386 on 11p12 within a TRAF6 intron (8.3×10−4 replication, 3.9×10−6 overall), and rs12746613 on 1q23.3 near FCGR2A (2.2×10−3 replication, 1.5×10−5 overall). These SNP associations likely represent true RA loci, but additional genotyping will be necessary for definitive confirmation.
Interestingly, many of the SNPs picked by GRAIL that validated in independent genotyping were not those with strongest evidence of association in the initial GWAS meta-analysis (see Figure 3A). That is, prioritization based purely on meta-analysis p-values might have missed many of these associations. For example, rs12746613 (FCGR2A) was ranked 163rd of 179 and rs540386 (RAG1/TRAF6) was ranked 110th. Of the 5 SNPs that we genotyped with the most significant GRAIL ptext scores, 3 replicated and 1 demonstrated nominal evidence of association; only rs2614394 (IRAK4) demonstrated no evidence of association.
Many of these seven alleles further implicate genomic regions already associated with autoimmune diseases (see Table 3). At this point none of these RA risk alleles correspond perfectly to any previously established autoimmune allele; but in some cases fine mapping of the region in multiple diseases could clarify the relationship between them. The rs12746613 FCGR2A SNP is 13 kb away from a missense SNP in FCGR2A that has been associated with systemic lupus erythematosus17,18; they are in the same LD block (r2=0.19, D’=1.0). The rs394581 SNP is located in the 5′ untranslated region of TAGAP and is 17 kb away from a SNP associated with celiac disease and with type I diabetes19,20; they are in partial LD (r2=0.32, D’=.73). The rs10919563 PTPRC SNP is 35 kb away from a rare (~1% allele frequency) non-synonymous SNP that alters splicing of PTPRC21; there have been inconsistent reports that it is associated with multiple sclerosis22-24. We also note that the rs7234029, a PTPN2 intronic SNP, is 41 kb away from a SNP associated with both type I diabetes and celiac disease20; they are in the same LD block (r2=0.14, D’=1.0). The rs548234 SNP is located 10 kb downstream from the PRDM1 transcript and is 133 kb away from a SNP previously associated with Crohn’s disease25. The rs11586238 SNP is 50 kb upstream of the CD2 start site, but is also near multiple other key immunological genes including CD58, and IGSF2. It is 159 kb away from a multiple sclerosis associated SNP within a CD58 intron26,27.
The rs1980422 SNP is located about 10 kb away from the 3′ UTR of CD28 and is 129 kb away from a known RA and type I diabetes risk allele in the CTLA4 region (rs3087243)11. There is minimal LD between these two SNPs (r2=0.04, D’=0.40); conditional analysis confirmed that the two SNPs independently confer RA risk (see Supplementary Table 5 online).
These SNP associations continue to clarify critical biological processes involved in RA pathogenesis, including T-cell activation, NF-κB signaling, and B-cell activation and differentiation. The CD2 protein is a co-stimulatory molecule on the surface of natural killer cells and T-cells; CD2 signaling is mediated by binding PTPRC directly28. Association to CD28, contributes additional evidence of the role of T-cell activation in disease pathogenesis. TRAF6 is involved in downstream NF-κB activation; it binds CD40 directly and is a key component of B-cell activation29. Our study has also implicated novel processes represented by PRDM1 (BLIMP-1), a transcription factor that regulates terminal differentiation of B-cells into immunoglobulin secreting plasma cells30. Functional studies and re-sequencing will be required to confirm that these genes are indeed the truly causal genes.
We examined all 7 replicated RA SNPs along with known RA risk alleles for epistatic interactions (see Supplementary Note online). Despite the functional relationships between these genes, we found no evidence of significant interactions.
Population stratification could result in spurious associations. However, we were careful for each collection either to use (1) epidemiologically matched samples or (2) ancestry informative markers to match cases and controls. We further note that our seven replicated SNP associations demonstrate consistent effects across all 14 collections without evidence of heterogeneity (p>0.05 by Breslow-Day test of Heterogeneity, see Table 2).
In this study we demonstrate the utility of functional information to prioritize SNPs for replication. We did not pre-define pathways, but rather we used GRAIL to look for genes that had relationships to other validated RA genes. We note that GRAIL is limited in its ability to identify disease genes in entirely novel pathways (i.e., pathways not suggested by the 16 previously known RA risk loci). Arguably, those loci could point to truly novel pathogenic mechanisms. Additionally, successful application of GRAIL is contingent on the scientific literature’s comprehensive description of relevant gene relationships. The general application of GRAIL to other diseases will depend critically on the completeness of the validated loci list and the documentation about relevant processes in the literature. Despite these limitations, our study has identified at least three novel RA risk loci, with strong evidence for additional risk loci.
GRAIL is a method that leverages statistical text mining principles to assess whether putative disease loci harbor genes with functional relationships to genes in other associated disease loci2. Two genes are considered similar if the words used to describe them in PubMed abstracts suggest similar functionality. The implementation of GRAIL used here leverages a text database of 250,000 abstracts published prior to December 2006.
To test the ability of GRAIL to distinguish RA risk loci from spurious associations we defined a set of true positive loci that were discovered since December 2006; these loci would not be described in the GRAIL text database. We also approximated a set of spurious associations by randomly selecting 10,000 SNPs from the Affymetrix 500K genotyping array. We tested both SNP sets for relationships to known associated RA loci with GRAIL. Validated SNPs were tested against the other 15 independent loci; spurious SNPs were tested against all 16 loci. The sensitivity was defined as the proportion of true positive associations that GRAIL assigned a ptext<0.01 score; the specificity was defined as the proportion of spurious associations that GRAIL assigned a ptext>0.01 score.
In order to identify SNPs for follow-up, we examined the results of a recently published meta-analysis of three GWAS studies (see Supplementary Table 3)1. We examined 336,721 SNPs outside the MHC region that passed strict quality control criteria. We identified those SNPs that were nominally associated with RA (p<0.001). We grouped SNPs into independent loci; two SNPs were placed in the same locus if there was evidence of LD (r2>0.1 in CEU HapMap). We removed all loci that overlapped with validated RA risk regions (see Table 1). We also removed loci with p<10−4 that were genotyped in most available patient collections and had failed to validate in a previous study1. From the remaining set of independent loci, we selected the single SNP that demonstrated the greatest evidence of association in the published meta-analysis.
We tested 179 candidate SNPs for relationships to genes within the 16 independent loci known to be associated with RA using GRAIL. SNPs that obtained compelling GRAIL scores (ptext<0.01) were selected for follow-up investigation. To assess the degree of enrichment among high scoring SNPs we sampled 100,000 random sets of 179 SNPs and tested these SNP sets with GRAIL. We calculate the proportion of sets with as many or more GRAIL hits to calculate a permutation based p-value. We note that the version of GRAIL that we used is a previous implementation that differs slightly from the published implementation2 - results are not substantially affected when we do the same experiment with the current version of GRAIL (see Supplementary Figure 2).
Patient collections are described in detail in Supplementary Table 3 and in the Supplementary Note. Each collection consisted only of individuals that were self-described white and of European descent, and all cases either met 1987 ACR classification criteria or were diagnosed by board certified rheumatologists. Informed consent was obtained from each patient, and the Institutional Review Board at each collecting site approved the study.
All cases were autoantibody positive (CCP and/or RF). For most of the collections, matched control samples were collected along with case samples as part of the same study. For some of the collections - where control samples were unavailable - we matched these case collections to shared controls. We used a total of eleven separate patient collections for replication genotyping: (1) CCP positive cases from the Brigham Rheumatoid Arthritis Sequential Study (BRASS)31 and controls from three separate studies on multiple sclerosis32, age-related macular degeneration33, and myocardial infarction34; (2) CCP positive cases from the Toronto area (CANADA)13 and controls recruited from the same site along with additional controls taken from a disease study of lung cancer35; (3) CCP positive cases and controls from Halifax and Toronto (CANADA-II)13; (4) CCP positive cases from Sweden and epidemiologically matched controls (EIRA-II)36; (5) CCP positive Dutch cases and controls collected from the greater Amsterdam region (GENRA)37,38; (6) North American RF positive cases and controls matched on gender, age, and grandparental country of origin from the Genomics Collaborative Initiative (GCI)4; (7) CCP or RF positive Dutch cases and controls from Leiden University Medical Center (LUMC) 39,40; (8) CCP positive cases drawn from North American clinics and controls from the New York Cancer Project (together this collection is called NARAC-II)13,36; (9) CCP positive cases drawn from North American clinics (NARAC-III)13 and publicly available controls taken from a Parkinson’s study41 and study 66 and 67 of the Illumina Genotype Control Database; (10) CCP or RF positive cases identified by chart review from the Nurses Health Study (NHS) and matched controls based on age, gender, menopausal status, and hormone use42; and (11) CCP or RF positive cases recruited at multiple sites in the United Kingdom by the United Kingdom Rheumatoid Arthritis Genetics (UKRAG) collaboration6. We used available SNP data from this and previous studies to identify genetically identical samples from the same country; we assumed these represented duplicated individuals and we removed them.
Detailed description of genotyping is provided in the Supplementary Note. All GWAS meta-analysis genotyping was previously described. We genotyped replication samples at the Broad Institute using a single Sequenom iPlex Pool (for EIRA-II and GENRA collections) and Affymetrix 6.0 (BRASS), the National Institutes of Health using a single Sequenom iPlex Pool (NARAC-II), Analytic Genetics Technology Centre in Toronto using a single Sequenom iPlex Pool (CANADA-II), Epidemiology Unit at The University of Manchester using a single Sequenom iPlex Pool (UKRAG), Celera using kinetic PCR43 (GCI and LUMC), at the Nurses Health Study in Boston using the BioTrove multiplex SNP genotyping assay (NHS), at the Feinstein Institute using the Illumina 317K array (NARAC-III); and at Illumina using the Illumina 370K array (CANADA). For NARAC-III we additionally obtained publicly available shared controls genotyped on a similar platform from two separate studies. In the cases where whole genome data were available we either extracted data for the 22 SNPs (BRASS) or used imputation to estimate genotypes for them (CANADA and NARAC-III).
For each collection we applied stringent quality control criteria. We required that each SNP pass the following criteria for each collection separately: (1) genotype missing rate < 10%, (2) minor allele frequency > 1%, and (3) Hardy-Weinberg equilibrium with p>10−3. We then excluded individuals with data missing for > 10% of SNPs passing quality control.
For each replication collection we corrected for possible population stratification by either (1) using only epidemiologically matched samples when cases and controls were drawn from the same population, or (2) matching at least one control for each case based on ancestry informative markers (see Supplementary Note for details). Since the cases in the NHS, GCI, LUMC, EIRA-II, CANADA-II, UKRAG, and GENRA collections were well matched to controls, we did not pursue further strategies to correct for population stratification. For the BRASS, NARAC-II, CANADA, and NARAC-III, we matched cases and controls with ancestry informative markers and placed them each into a single stratum. For the BRASS cases and shared controls, GWAS data on Affymetrix 6.0 (unpublished data) was available; we used 681,637 SNPs passing strict quality control as ancestry informative markers. For NARAC-II cases and NYCP shared controls, case and controls were matched using genotype data on 760 ancestry informative markers. For the NARAC-III cases and shared controls, we used available Illumina 317K GWAS data for 269,771 SNPs passing stringent quality control criteria. For the CANADA cases and controls, we used available Illumina 317K GWAS data for 269,771 SNPs passing stringent quality control criteria. For each case-control collection, we used these SNPs to define the top 10 principal components and to remove genetically-distinct outliers (sigma threshold = 6 with five iterations) with the software program EIGENSTRAT44. We eliminated vectors that correlated with known structural variants on chromosomes 8 and 17, demonstrated minimal variation, or did not stratify cases and controls. After mapping cases and controls in the space of eigenvectors, we matched cases to controls that were nearest in Euclidean distance as described elsewhere1.
For each SNP we conducted three statistical tests. First, we conducted a one-sided CMH statistical test across eleven strata to assess if RA association was reproducible in the replication collections in the same direction as the GWAS meta-analysis. We set our significance threshold, after correcting for 22 hypothesis tests, to be p<0.0023 (=0.05/22). Second, we conducted a 573 strata joint analysis across all meta-analysis strata and substrata and replication strata; the eleven replication collections were each placed into their own strata, while the meta-analysis samples were partitioned into 562 strata to be consistent with the approach taken in the original analysis to correct for stratification1,36. Third, we calculated a Breslow-Day test of heterogeneity of odds ratios. We performed all analyses in MATLAB.
SR is supported by an NIH Career Development Award (1K08AR055688-01A1) and an American College of Rheumatology Bridge Grant. RMP is supported by a K08 grant from the NIH (AI55314-3), a private donation from the Fox Trot Fund, the William Randolph Hearst Fund of Harvard University, the American College of Rheumatology ‘Within Our Reach’ campaign, and holds a Career Award for Medical Scientists from the Burroughs Wellcome Fund. MJD is supported by NIH grants through the U01 (HG004171, DK62432) and R01 (DK083756-1, DK64869) mechanisms. The Broad Institute Center for Genotyping and Analysis is supported by grant U54 RR020278 from the National Center for Research Resources. The BRASS Registry is supported by a grant from Millennium Pharmaceuticals and Biogen-Idec. The NARAC is supported by the NIH (NO1-AR-2-2263 and RO1 AR44422). This research was also supported in part by the Intramural Research Program of the National Institute of Arthritis, Musculoskeletal and Skin Diseases of the National Institutes of Health. This research was also supported in part by grants to KAS from the Canadian Institutes for Health Research (MOP79321 and IIN - 84042) and the Ontario Research Fund (RE01061) and by a Canada Research Chair. We acknowledge the help of C.Ellen van der Schoot for healthy control samples for AMC/UVA and the help of Ben A.C. Dijkmans, Dirkjan van Schaardenburg, A. Salvador Peña, Paul L. Klarenbeek, Zhuoli Zhang, Mike T Nurmohammed, Willem F Lems, Rob R.J. van de Stadt, Wouter H. Bos, Jenny Ursum, Margret G.M. Bartelds, Daniëlle M. Gerlag, Marleen G.H. van der Sande, Carla A. Wijbrandts, and Marieke M.J. Herenius in gathering GENRA patient samples and data. We thank the Myocardial Infarction Genetics Consortium (MIGen) study for the use of genotype data from their healthy controls in our study. The MIGen study was funded by the U.S. National Institutes of Health and National Heart, Lung, and Blood Institute’s STAMPEED genomics research program R01HL087676 and a grant from the National Center for Research Resources. We thank Johanna Seddon Progression of AMD Study, AMD Registry Study, Family Study of AMD, The US Twin Study of AMD, and the Age-Related Eye Disease Study (AREDS) for use of genotype data from their healthy controls in our study. We thank David Hafler and the Multiple Sclerosis collaborative for use of genotype data from their healthy controls recruited at Brigham and Women’s Hospital.