Following the identification of a disease-associated region by GWAS, comprehensive study of sequence variation in the region is required to identify the full set of variants that might explain the association signal. Since GWAS arrays incompletely capture DNA variation in each region, it has been hypothesized that causal variants partially captured by linkage disequilibrium (LD) – due to location near recombination hotspots or lower minor allele frequency – might, if directly tested, display stronger association to phenotype than the tag SNPs used in GWAS. In particular, because HapMap and GWAS arrays contain primarily variants with minor allele frequency (MAF) >5%, first generation GWAS studies failed to test polymorphisms of somewhat lower frequency that might have larger effects on disease risk. Finally, even in regions where the true association signal is well captured by LD to array SNPs, enumeration of all associated variants is a necessary prerequisite to functional experiments that will identify causal mutation(s). Thus, an important next step following GWAS is to assemble a more complete catalog of variation present in an associated region, and to test it for association to the phenotype of interest.

With the advent of next generation sequencing and the emergence of data from the 1000 Genomes (1000G) Project, investigators must choose between (or combine) multiple strategies for creating and testing a reference panel of polymorphic sites. We re-sequenced ~240kb on chromosome 9p21 (chr9:21936711-22176221, hg18) spanning the T2D and MI associations in 47 unrelated individuals of European ancestry from the HapMap CEU population¹⁷ as part of a sequencing project spanning six T2D-associated regions (Supplementary Table 1). Sequencing was performed at the Broad Institute on Illumina Genome Analyzers (Supplementary Note, all data available in the NCBI Short Read Archive). An analytical framework (Supplementary Note, Supplementary Table 2, Supplementary Figs. 1-5), since extended and incorporated in the Genome Analysis Tool Kit^18,19, was developed and includes methods to empirically recalibrate Illumina base quality scores, a Bayesian framework to call SNPs, local re-alignment to identify insertions/deletions (and remove clusters of false positive SNPs), and filters to remove false positive SNP calls based on discrepancy between forward and reverse strands.

This targeted sequencing identified 635 high-confidence SNPs on chromosome 9p21 (4,463 across the six regions) (Supplementary Table 3, Supplementary Fig. 6, SNPs available in dbSNP). We evaluated sensitivity against HapMap II¹⁷ and the high coverage Pilot 2 data from the 1000 Genomes Project²⁰ (Supplementary Note): at sites in overlapping samples with 10x or greater read coverage (70% of the region), sensitivity was 99% for HapMap variants and 97% for variants found in 1000G Pilot 2 (Supplementary Fig. 7a-c). To evaluate specificity, we genotyped 257 sites found on chromosome 9p21 but not previously genotyped in HapMap (Supplementary Fig. 7d, e). Overall, 96% of variants seen more than once in sequencing validated in the genotyping data (Supplementary Table 4).

We compared these variants to those discovered in the low-coverage Pilot 1 of the 1000 Genomes Project ²⁰, limiting comparison to 32 CEU individuals studied in both projects. Across the six regions, both projects identified similar numbers of variants: 3,897 SNPs in the high coverage targeted sequencing as compared to 4,043 in 1000G Pilot 1. However, the variants found were in fact only partially overlapping. Of variants seen in the high coverage targeted sequencing, 22% were missed by 1000G Pilot 1 (Fig. 1), nearly all of which were rare: 72% of these sites were singletons and 12% were seen twice (Fig. 1, Table 1). (Pilot 1 successfully identified 97% of SNPs seen more than 5 times in high coverage sequencing (Table 1)). Of variants identified in Pilot 1 but not in targeted sequencing (n=998), nearly all were sites at which target capture failed to achieve high coverage: 65% of these sites had zero coverage. Thus, targeted capture and low-pass whole genome had distinct and non-overlapping failure modes.

Open in a separate window

Figure 1

Comparison of targeted sequencing to 1000G Pilot 1 Data

Variant calls were made in all six regions of T2D association in the 32 individuals who were sequenced as part of both this targeted, high coverage sequencing effort (total 47 CEU HapMap individuals) as well as 1000G Pilot 1 (total 60 CEU HapMap individuals).

We evaluated methods for testing these variants for association to disease via linkage disequilibrium and haplotype-based imputation. First, we genotyped SNPs found in targeted re-sequencing on chromosome 9p21 in 168 individuals (56 parent offspring trios) from the HapMap extended CEU population²¹ (Supplementary Note). We used MACH^22,23 to impute variants from this reference panel into 1,000 T2D patients and 1,048 controls from the Diabetes Genetics Initiative (DGI) cohort¹ and 1,274 MI cases and 1,407 controls from the Myocardial Infarction Genetics (MIGen) Consortium cohort⁶, each previously genotyped on Affymetrix GWAS arrays (Supplementary Note).

We compared the results of imputation with this augmented reference panel (n=464 variants, Supplementary Table 5) to those obtained when imputing from HapMap II alone (n=238 variants). The addition of genotype data for a more complete collection of common variants provided imputation data for a much larger number of SNPs than was possible with HapMap II, which contains only 50-60% of common variants (Fig. 2a, b and Supplementary Fig. 8a, b). However, even with the augmented reference panel, the tag SNP density characteristic of the first generation GWAS arrays on which our disease samples were typed allowed only 80% of common (MAF > 5%) variants to be captured (either directly typed or imputed with a MACH-estimated r² ≥ 0.8). Moreover, only a small fraction of intermediate frequency variation (MAF 2-5%) was imputed with an estimated r² above this stringent threshold (Fig. 2c, d and Supplementary Fig. 8c, d).

Open in a separate window

Figure 2

Fraction of variation on chromosome 9p21 captured in T2D disease cohort by different imputation scenarios

MACH imputation quality estimates (a, c, e) and overall fraction of variation captured in T2D samples (b, d, f) for different imputation scenarios. (a, c, e) The MACH-estimated r² for each SNP is plotted as a function of genomic position. SNPs not observed in the reference panel are assigned an r² of zero. Recombination rate (estimated from HapMap) is plotted to reflect local LD structure. Gene annotations were taken from the University of California-Santa Cruz Genome Browser. (b, d, f) The fraction of variants captured in T2D samples is shown as a function of MAF and MACH-estimated r². Imputation scenarios are: (a, b) Imputing from HapMap II (n=238 SNPs in 60 individuals) into the SNPs genotyped on the Affymetrix 500K array; (c, d) Imputing from 112 individuals genotyped at HapMap II sites and validated sequencing sites (total n=464 SNPs) into the SNPs genotyped on the Affymetrix 500K array; (e, f) Imputing from the same reference panel as c, d into SNPs genotyped on the Affymetrix 500K array plus additional tag SNPs genotyped in the T2D cohort (genotyped marker density in T2D samples ~1 SNP/1.5kb).

To evaluate the impact of tag SNP density on imputation performance, we increased the number of tags across the region to approximately 1 SNP per 1.5kb (the previous density was ~1SNP/5kb in T2D samples and ~1SNP/3kb in MI samples) in the T2D and MI cohorts (Supplementary Note). With this increased density of tag SNPs, nearly all common variants (~98%) were captured with r² ≥ 0.8 in disease samples. Moreover, performance for intermediate frequency variants was dramatically improved, rising from 2% to 75% with r² ≥ 0.8 (Fig. 2e, f and Supplementary Fig. 8e, f). This result was not specific to the Affymetrix GWAS arrays, as we observed a similar improvement in imputation ability upon addition of tag SNPs using multiple other GWAS arrays (Supplementary Fig. 9).

We next compared different reference panels, imputing in each case into disease samples with the higher tag SNP density. The reference panels were: (a) the genotyped reference panel of 168 individuals above (112 unrelated individuals), (b) the targeted sequencing data (47 individuals, without genotyping and expansion into a larger sample set), and (c) 1000 Genomes Pilot 1 (55 individuals). We considered both the fraction of variants in each reference panel successfully imputed (which is related to the quality and completeness of SNP genotypes and to the size of the reference panel) and the fraction of all variation captured (which, in addition, depends on the proportion of all known SNPs represented in the reference panel).

The union of the three reference panels contained 582 variants (Fig. 3a). Each panel was partially incomplete, due to genotyping assay failure in the genotyped panel (14% of SNPs missing), sample size and low coverage in 1000 genomes (16% of SNPs missing), and sample size and gaps in coverage in the targeted sequencing (19% of SNPs missing). For common variants, there is little difference in bulk performance between the reference panels. Considering only SNPs contained in each reference panel (Fig. 3b) the genotyped panel has the highest proportion of variants imputed well. However, when all variation is considered (Fig. 3c), a lower proportion of common variation is captured by imputing from the genotyped reference panel, owing to the fact that some SNPs were missing in this panel because they failed assay design or quality control. Notably, the 1000G (freely available) and sequencing (costly) strategies performed equivalently for these common variants.

Open in a separate window

Figure 3

Comparison of imputation from a genotyped reference panel, directly from high coverage re-sequencing data, and directly from 1000G Pilot 1 data

(a) Variants present in the three reference panels and their overlap. The 67 variants present in the genotyped reference panel but not in the high coverage sequencing reference panel (denoted by asterisk) were called in high coverage sequencing as singletons and so were excluded from the sequencing reference panel. 40% of these variants are not singletons in the larger genotyped reference panel. (b) The fraction of sites within each reference panel captured with a MACH-estimated r² of at least 0.8. (c) The overall fraction of known variants captured with a MACH-estimated r² of at least 0.8 by imputation from each reference panel.

For intermediate frequency variants, there are more pronounced differences between the panels (Fig. 3b, c). These variants were best imputed from the genotyped reference panel (Fig. 3b), which was the largest and also contained trio information. This was true even when all variation was considered (Fig. 3c), suggesting that the improved imputation quality from genotype data and increased sample size offset the loss of variants in this panel due to genotyping failure. Comparing the high coverage re-sequencing and 1000G reference panels, lower frequency variants were better imputed from the high coverage re-sequencing data both when considering only the SNPs within each reference panel (Fig. 3b) and when considering the overall proportion of low frequency variants captured by imputation from each reference panel (Fig. 3c). This is consistent with the low coverage 1000G pilot 1 data being less complete and accurate for lower frequency variants²⁰.

We tested variants for association to T2D and MI using imputation from all three reference panels to maximize the number of variants captured (Supplementary Note). Overall, we have captured 461 of the 582 polymorphic variants observed across all three reference panels in our T2D and MI samples with a MACH-estimated r² of at least 0.8: this represents ~ 92% of all known common variants and ~52% of intermediate frequency variants (at a MACH-estimated r² of 0.5, these figures are 98% of common variants and 83% of intermediate frequency variants). In comparison, only 176 of the 582 variants were previously captured by imputation from HapMap. Despite testing many additional SNPs in partial LD with the index GWAS hits and at allele frequencies not well captured by first generation GWAS arrays and HapMap, we found no example of a SNP with stronger association to T2D or MI than the initial GWAS signals.

However, we did identify multiple additional variants in strong LD with the GWAS hits that might underlie each association. We observed the three-tiered haplotypic association to T2D first reported by the Wellcome Trust Case Control Consortium with protective, risk, and neutral haplotypes (Table 2). The protective alleles of the GWAS SNP (rs10811661) and nine other SNPs in strong LD with this variant tag the protective haplotype (Fig. 4a, Supplementary Table 6). Interestingly, no single SNP yet identified marks the risk haplotype. Association analyses for MI identified 7 SNPs in LD with each other and with equivalent evidence for association (P < 10⁻⁴) as well as 54 additional SNPs with only slightly less evidence for association (P < 10⁻³) (Fig. 4b, Supplementary Table 6). Knockout of the MI-associated region in mouse alters regulation of CDKN2A and CDKN2B²⁴, and two of the associated SNPs have recently been shown to disrupt a STAT1 binding site²⁵. Interestingly, in addition to the SNPs disrupting the STAT1 site, there are other variants with equivalent MI association and with putative functional annotations, including variants overlapping exons of the non-coding transcript CDKN2BAS, highly conserved regions, and predicted, conserved transcription factor binding sites (Supplementary Table 6).

Open in a separate window

Figure 4

Association results for T2D and MI on chromosome 9p21

Regional plots showing association signal for (a) T2D and (b) MI. The signal for each SNP (represented as −log₁₀ p-value) is plotted as a function of genomic position. The size of the diamond for each SNP represents the LD (measured as r²) between that SNP and the original GWAS SNP (rs10811661 for T2D and rs4977574 for MI). Recombination rate (estimated from HapMap) is plotted to reflect the local LD structure in the region. Gene annotations were taken from the University of California-Santa Cruz Genome Browser.

This study is limited by the investigation of a single region (albeit one with at least eight different disease associations), by the early nature of the sequencing data analyzed, by the small number of samples sequenced in SNP discovery, and by the sample size of our disease cohorts. Nonetheless, the observations on the strengths and weaknesses of different methods for fine mapping GWAS signals are likely general: targeted high coverage sequencing provides high sensitivity for lower frequency variants, but has gaps in coverage; the 1000G Pilot 1 resource offers more even coverage at lower depth, currently sufficient for capture of most common variation; creating a genotyped reference panel improves accuracy and sample size, but is limited by assay conversion failures; tag SNP density characteristic of first generation GWAS is inadequate to maximally extract information with current imputation algorithms. To some extent, these limitations are transient: the growing 1000 Genomes Project resource is sequencing over 2,000 diverse samples with both low-pass whole genome and high coverage targeted exon approaches, increasing the accuracy and completeness of the public reference panel. However, our results suggest that fully exploiting this resource for imputation may require increasing tag density in GWAS disease samples and / or improved algorithms for imputation.

Finally, our study did not find evidence for stronger association at 9p21 to SNPs in moderate LD with the initial tags. While the maximum achievable association signal for lower frequency variants was limited by our sample size, we did not observe lower frequency variants with effect sizes that could individually explain the common variant associations. We do, however, identify additional common variants in LD with the GWAS hits that might underlie each association. Enumeration of all variants on 9p21 that might explain each association signal will be needed as a foundation for systematic functional studies that aim to understand how different non-coding variants in this single genomic interval can lead to such varied and clinically significant phenotypic associations.