|Home | About | Journals | Submit | Contact Us | Français|
Genome-wide association studies (GWAS) conducted using commercial single nucleotide polymorphisms (SNP) arrays have proven to be a powerful tool for the detection of common disease susceptibility variants. However, their utility for the detection of lower frequency variants is yet to be practically investigated. Here we describe the application of a rare variant collapsing method to a large genome-wide SNP dataset, the Wellcome Trust Case Control Consortium rheumatoid arthritis (RA) GWAS. We partitioned the data into gene-centric bins and collapsed genotypes of low frequency variants (defined here as MAF ≤0.05) into a single count coupled with univariate analysis. We then prioritised gene regions for further investigation in an independent cohort of 3,355 cases and 2,427 controls based on rare variant signal p value and prior evidence to support involvement in RA. A total of 14,536 gene bins were investigated in the primary analysis and signals mapping to the TNFAIP3 and chr17q24 loci were selected for further investigation. We detected replicating association to low frequency variants in the TNFAIP3 gene (combined p = 6.6 × 10−6). Even though rare variants are not well-represented and can be difficult to genotype in GWAS, our study supports the application of low frequency variant collapsing methods to genome-wide SNP datasets as a means of exploiting data that are routinely ignored.
The online version of this article (doi:10.1007/s00439-010-0889-1) contains supplementary material, which is available to authorized users.
Genetic epidemiology of common disease has been transformed by the recent advent of the genome-wide association study (GWAS) that allows the unbiased testing of hundreds of thousands of single nucleotide polymorphisms (SNPs) in large sample collections. This approach has greatly accelerated the identification of genetic markers associated with disease susceptibility and continuous traits. The Catalogue of Published Genome-Wide Association Studies contains details from 378 published GWAS reporting 1,700 SNP associations to diseases or quantitative traits (with p < 10−5) to date (Hindorff et al. 2009).
While the GWAS strategy has undoubtedly been successful in identifying a multitude of disease-associated makers, these have been shown to explain very little of the estimated heritable component of complex disease (Maher 2008). For example, all of the currently known rheumatoid arthritis (RA) risk alleles, including the major histocompatibility complex (MHC) locus, have been estimated to explain less than half of the total genetic component (Raychaudhuri et al. 2008), indicating that a substantial proportion of susceptibility variants remain undiscovered.
The deficit between estimated heritability of a particular trait and its known genetic variation may be attributed to a number of issues. Firstly, the SNPs identified in the primary scans are likely to serve only as proxies for the true causal variant, therefore estimates of their contribution to genetic risk may be greatly underestimated (McCarthy and Hirschhorn 2008). Secondly, further variants at an associated locus may confer additional effects independent of the primary signal (Plenge et al. 2007). In addition, structural variation, such as insertions, deletions, and duplications, is now recognised as a well-established feature of the genome (Iafrate et al. 2004; Redon et al. 2006; Sebat et al. 2004), with a number of emerging studies reporting associations to disease susceptibility (Fanciulli et al. 2007; Fellermann et al. 2006; Gonzalez et al. 2005; Hollox et al. 2008; McKinney et al. 2008). Importantly, a significant proportion of inherited susceptibility to common disease may be attributable to the effects of multiple, low frequency variants with moderate effects in a number of genes (Bodmer and Bonilla 2008). The majority of variants in the human genome are predicted to have minor allele frequency (MAF) below 0.05 and are correlated with increased functional significance (Gorlov et al. 2008). Evidence to support the involvement of rare variants in the susceptibility to common disease comes from both theoretical (Pritchard 2001; Pritchard and Cox 2002) and empirical data from studies such as colorectal cancer (Fearnhead et al. 2004), HDL cholesterol levels (Cohen et al. 2004), hypertension (Ji et al. 2008), obesity (Ahituv et al. 2007) and type 1 diabetes (Nejentsev et al. 2009).
The design of currently available commercial SNP platforms has been heavily biased towards common SNPs. GWAS rely on indirect association methods exploiting local linkage disequilibrium (LD) patterns. Studies performed on large empirical and simulated datasets have shown that LD-tagging based on HapMap data is robust for common SNPs, but that performance declines rapidly for lower frequency variants when common tags are used (Ahmadi et al. 2005; Zeggini et al. 2005). In addition, rare variants are notoriously difficult to genotype accurately, and single-point association analysis of rare variants has inherently low power given currently available sample sizes. Indeed, the majority of GWAS exclude low frequency SNPs from further analysis, thus ignoring any information afforded by the rare variants included in the arrays. The issue of power to detect association is being addressed by the development of statistical methods for the combined analysis of rare variants (Li and Leal 2008; Morris and Zeggini 2010).
Here we demonstrate, with a proof of principle experiment, that the direct analysis of rare variants (defined here as MAF ≤0.05) present on commercially available arrays can indeed lead to the identification of real complex disease susceptibility signals. We have analysed all rare variants on the Affymetrix 500K array for the RA branch of the Wellcome Trust Case Control Consortium (2007) and have followed up four signals with strong priors for association with disease. Following extensive checks on genotype quality and after excluding poorly clustering variants, we find replicating association of rare variation at the TNFAIP3 locus with RA susceptibility.
This study was approved by the North West Multicentre Research Ethics Committee (MREC 99/8/84), and all subjects provided informed consent.
We used data from the RA branch of the WTCCC study (The Wellcome Trust Case Control Consortium 2007), genotyped on the Affymetrix GeneChip Human Mapping 500K Array Set. The dataset consisted of 1,860 RA cases and 2,938 controls, all from the UK, after the exclusion of samples failing quality control (QC) (The Wellcome Trust Case Control Consortium 2007). We excluded from further analysis SNPs failing QC thresholds as described in the WTCCC study (The Wellcome Trust Case Control Consortium 2007). Briefly, SNPs were excluded based on deviation from Hardy–Weinberg equilibrium (p < 5.7 × 10−7) or if the study-wise missing data proportion was >0.05. Additionally, SNPs with a MAF <0.05 were excluded if the study-wise missing data proportion exceeded 0.01. Genotype cluster plot inspection was not possible at the genome-wide scale, so we inspected the clustering properties of SNPs within signals of interest after the rare variant analysis had been carried out.
In order to maximise power to detect association, we followed a “super-locus” approach to analysing rare variation in the GWAS RA data (Li and Leal 2008; Morris and Zeggini 2010). SNPs with a study-wise MAF <0.05 were included in the analysis. We defined genic regions genome-wide based on the coordinates of known genes and (arbitrarily) included 50 kilobases (kb) flanking either side of each gene’s transcriptional start and stop site in an attempt to include SNPs affecting localised regulatory elements. For each such region, we collapsed the rare variant allele counts by assigning each individual a carrier/non-carrier status based on the presence or absence of rare variant minor alleles (Supplementary Figure 1) (Li and Leal 2008; Morris and Zeggini 2010). This classification system is not affected by LD, as the number of rare alleles carried by each individual (the rare allele load) was not taken into account. This allowed the construction of a simple 2 × 2 contingency table, which we used to test differences in rare variant minor allele carriage between cases and controls using the Chi-squared test (or Fisher’s exact test where necessary). For regions attaining p < 10−4, we also permuted case/control status 100,000 times in order to assess significance by means of empirical p values. We used CCRaVAT software to carry out these analyses (http://www.sanger.ac.uk/resources/software/ccravat-qutie/). For the purposes of this study, variants with a MAF ≤0.05 were considered low frequency/rare and included in the analysis. Although this value is arbitrary, it represents a typical threshold at which SNPs are excluded from further consideration in many GWAS.
Signals were selected for replication based on the primary rare variant association scan p value and the existence of prior evidence to suggest involvement in the susceptibility to RA. Bins containing only a single rare variant were excluded from further investigation, as were bins mapping to the extended MHC region. For the purposes of this study, the extended MHC was considered to be a 7.6 Mb region consisting of the extended class I, classical class I, extended class II, classical class II, and classical class III subregions. The region is defined by the co-ordinates Chr6:25,809,997-33,486,772 (NCBI build 36). As low frequency/rare variants are notoriously difficult to accurately genotype, we examined SNP clustering properties for all variants contributing to the prioritised signals. To investigate the influence of poor quality genotyping on the performance of the super-locus method we excluded any SNPs considered to be of low quality, based on clustering, from the gene region under study and re-analysed.
DNA samples for 3,838 RA cases and 2,719 healthy unrelated controls were made available from the UK Rheumatoid Arthritis Genetics (UKRAG) consortium, a collaboration of six UK rheumatology research groups.
Low frequency/rare variants contributing to prioritised association signals in the primary scan were selected for follow-up. Genotyping was performed using the Sequenom MassARRAY platform with 10 ng of genomic DNA in concordance with the manufacturer’s specifications for the iPLEX protocol. SNPs rs5029939 and rs7749323 had already been typed in a subset of the follow-up patient samples as part of a related study (Orozco et al. 2009). All Sequenom genotype cluster plots were manually evaluated to determine satisfactory resolution of clusters. Samples with a call rate <0.95 and SNPs with a call rate <0.90 were excluded prior to analysis. Allele counts for all SNPs mapping to a gene were combined and then compared between cases and controls using a Fisher’s exact test.
We combined individual summary statistics across the WTCCC RA and UK replication studies using a fixed effects inverse-variance meta-analysis approach for each genic region. We also investigated the presence of significant heterogeneity across studies, by means of the I2 and Q statistics.
Following QC procedures, the GWAS consisted of 1,860 RA cases and 2,938 healthy controls genotyped for 459,446 SNPs. This dataset included 40,482 SNPs with MAF ≤0.05 suitable for inclusion in the current analysis. Application of the super-locus method yielded a total of 22,344 genomic bins, of which 14,536 contained at least two rare variants and mapped outside of the defined MHC region (λGC = 1.09, Supplementary Figure 2). We identified 25 gene regions with a p < 10−4. For the purposes of this study, we prioritised four signals for follow-up genotyping that emerged from two distinct chromosomal regions. This selection was based on their rare variant analysis p value and prior evidence to support disease association. Firstly, we selected SNPs from three gene regions mapping to chromosome 17q24; ACE, CYB561, and PRKCA (p = 6.0 × 10−5, p = 2.7 × 10−5 and p = 0.0012 respectively). This region is implicated in susceptibility to RA by previous evidence from human and animal linkage studies and is of particular interest to our research group (Backdahl et al. 2008; Barton et al. 2001; Jawaheer et al. 2001). To date, no common susceptibility variants that could explain the observed linkage signal have been discovered in this region, and this could be attributable to the presence of associated rare variants. The second candidate locus was the TNFAIP3 gene (p = 4.2 × 10−5), a known RA-associated locus with evidence for allelic heterogeneity (Orozco et al. 2009; Plenge et al. 2007; Thomson et al. 2007). These four bins comprise a total of 27 unique SNPs (due to overlap the ACE and CYB561 gene regions share two SNPs). Individual SNP genotype data and association results are reported in Supplementary Table 1. For the purposes of evaluating data quality we examined the clustering properties of rare variants residing in the associated gene regions and identified six SNPs with unsatisfactory clustering. Re-analysis of the data after excluding the badly clustering variants did not alter the results for the TNFAIP3 and PRKCA loci qualitatively, but did affect significance at the other two genes (Table 1; Supplementary Table 2).
We attempted to genotype all 27 SNPs from the four bins, regardless of clustering results, in an independent dataset consisting of 3,355 cases and 2,427 controls after QC filtering. Three SNPs failed assay multiplexing and 24 SNPs were successfully genotyped. SNPs were partitioned into appropriate gene regions and proportions of carriers and non-carriers were compared between the case and control groups. Evidence to support rare variant association with RA in the TNFAIP3 gene was further substantiated (p = 0.03) (Table 1). No evidence was found to support association in the three genes mapping to the chr17 candidate region.
Meta-analysis across the discovery and replication datasets revealed stronger evidence for association of low frequency variants in TNFAIP3 with RA (combined p = 6.6 × 10−6; OR 1.27[1.14:1.40]) (Table 1). The direction of association indicated an over-representation of rare variant minor alleles in cases compared to controls.
To date GWAS have concentrated on the investigation of common variation therefore discarding a significant proportion of data already collected, since a common QC criterion is the exclusion of SNPs with low MAF. The evidence presented here highlights the potential importance of investigating low frequency variant data from commercial SNP arrays. It should be noted that the primary purpose of these analyses is signal detection and not accurate effect estimation. The application of a rare variant analysis method (Li and Leal 2008; Morris and Zeggini 2010), which accumulates information from multiple rare SNP alleles within a genic region, identified an association between RA and the TNFAIP3 gene.
Multiple independent RA associations have previously been reported at the TNFAIP3 locus. The most convincing evidence for association to rare variants in this region comes from the meta-analysis described here (p = 6.6 × 10−6). Association to the 6q23 region was originally identified for common variants by two independent GWAS (Plenge et al. 2007; The Wellcome Trust Case Control Consortium 2007). The WTCCC RA GWAS identified a common RA risk variant >180 kb away from the TNFAIP3 gene, which was replicated in a UK study (Thomson et al. 2007). This finding prompted a fine mapping experiment covering the intergenic region containing the index association and the TNFAIP3 gene across ~7,000 individuals, identifying an independent association at rs5029937 (MAF ~0.04) in intron 2 of TNFAIP3 (Orozco et al. 2009). The evidence for multiple independent effects in the TNFAIP3 gene region make this a compelling candidate for re-sequencing.
Application of the super-locus approach using the original WTCCC data alone (2,000 cases, 3,000 controls) prior to the fine-mapping project would have directed attention to this signal obviating the need for time-consuming and costly region-wide fine-mapping. LD between the low frequency variants in the TNFAIP3 gene region reported here and the intergenic SNPs, previously shown to be associated with RA, is minimal (r2 < 0.04). However, there is strong LD between two SNPs from the super-locus analysed here and rs5029937 (rs7749323 r2 = 0.91 and rs5029939 r2 = 0.97). Our results highlight the potential strength of analysing lower frequency/rare variants in a complementary approach to the primary GWAS of common variants.
The study of low frequency/rare variants from commercial SNP arrays is limited by a number of factors. Firstly, the arrays contain only a small proportion of the low MAF variants that actually exist; this is due to the array design being motivated by the detection of common susceptibility variants. The gold standard strategy for investigating lower frequency/rare variants will be the use of sequencing data, capturing a much higher proportion of genetic variation (Morris and Zeggini 2010). However, in spite of the rapid advances in sequencing technology, this approach is often limited by economic constraints when considering large sample sizes. The method applied here utilises existing data that are generally discarded in GWAS. It therefore increases the amount of available data that are analysed and comes at no additional cost to the researcher. Collapsing methods are also very flexible in their application. In the study presented here we have focused on gene-centric bins, but this could be expanded to include bins based on coordinates of known conserved non-coding regions or functional pathways as the functional unit for analysis.
Secondly, the genotype quality of rare variants typed on GWAS platforms tends to be low, mainly driven by poor automated clustering and genotype calling. This limitation affects both single-point and collapsing methods and necessitates stringent and exhaustive inspection of cluster plots in association signals of interest. The potential for confounding is demonstrated in our study by the effect of excluding poorly clustering SNPs in the ACE and CYB561 gene regions. These two bins were no longer statistically significant upon removal of poor quality SNPs, thus highlighting the potential for inflated false positive rates in the absence of stringent quality assurance when dealing with low frequency SNP data. In addition, it is highly likely that bins containing only a few SNPs have an increased sensitivity to such artefacts. Genotype calling algorithms fine-tuned for rare variant clustering have recently been developed and hold the promise of helping to overcome many of these issues.
Thirdly, collapsing methods can be limited by the properties of the component variants. For example, they have been shown to be sensitive to the inclusion of neutral alleles, whose presence reduces power. Ideally, each of the bins would only contain putatively functional variants, which could be achieved by the use of bioinformatics tools. However, the accuracy of such tools is currently limited and the reduction of power is greater when true functional SNPs are excluded compared to the inclusion of functionally neutral SNPs (Li and Leal 2008). Additionally, it is conceivable that bins could contain variants with diverse effects, where the presence of both protective and risk alleles would confound the analysis. Rare variant collapsing methods are being further developed and extended to account for direction of effect, probability that a variant be functional, and genotype or sequence call uncertainty.
Finally, the phenotype under investigation could strongly dictate the success of identifying novel, rare disease-related variants. It could be argued that late-onset diseases, such as RA, are less likely to be influenced by a contribution from rare variants. However, the underlying genetic architecture of the disease is yet to be fully elucidated; collapsing methods represent one strategy that can help address this question.
Despite these limitations, our results illustrate that collapsing methods represent a strategy with the potential for identifying novel variants at no additional cost to the researcher, thus maximising return from costly primary GWAS. Applying this approach to existing SNP array data would have expedited the identification of the TNFAIP3 intron 2 SNP association with RA. We advocate the application of this method in conjunction with primary GWAS for complex diseases, but also suggest extreme caution with the interpretation of promising signals, which need to be scrutinised for genotype quality. Additionally, collapsing methods hold great promise for the future with the application to sequencing data.
Below is the link to the electronic supplementary material.
The authors are grateful to Edward Flynn and Paul Martin for Sequenom SNP genotyping. We are also grateful for the support of the Manchester Biomedical Research Centre and the Manchester Academy of Health Sciences.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.
United Kingdom Rheumatoid Arthritis Genetics (UKRAG). School of Medicine and Biomedical Sciences, Sheffield University, Sheffield S10 2JF (Dr Anthony G Wilson); NIHR-Leeds Musculoskeletal Biomedical Research Unit, Leeds Institute of Molecular Medicine, University of Leeds, UK (Prof. Ann W Morgan, Prof. Paul Emery); Clinical and Academic Rheumatology, Kings College Hospital NHS Foundation Trust, Denmark Hill, London SE5 9RS (Dr. Sophia Steer); Musculoskeletal and Genetics Section, Division of Applied Medicine, University of Aberdeen, UK, AB25 2ZD (Dr Lynne J Hocking, Dr David M Reid); University of Oxford Institute of Musculoskeletal Sciences, Botnar Research Centre, Oxford OX3 7LD, UK (Dr Pille Harrison, Professor Paul Wordsworth).
The members of the UKRAG are given in Appendix.