The Metabochip, a Custom Genotyping Array for Genetic Studies of Metabolic, Cardiovascular, and Anthropometric Traits

Genome-wide association studies have identified hundreds of loci for type 2 diabetes, coronary artery disease and myocardial infarction, as well as for related traits such as body mass index, glucose and insulin levels, lipid levels, and blood pressure. These studies also have pointed to thousands of loci with promising but not yet compelling association evidence. To establish association at additional loci and to characterize the genome-wide significant loci by fine-mapping, we designed the ''Metabochip,'' a custom genotyping array that assays nearly 200,000 SNP markers. Here, we describe the Metabochip and its component SNP sets, evaluate its performance in capturing variation across the allele-frequency spectrum, describe solutions to methodological challenges commonly encountered in its analysis, and evaluate its performance as a platform for genotype imputation. The metabochip achieves dramatic cost efficiencies compared to designing single-trait follow-up reagents, and provides the opportunity to compare results across a range of related traits. The metabochip and similar custom genotyping arrays offer a powerful and cost-effective approach to follow-up large-scale genotyping and sequencing studies and advance our understanding of the genetic basis of complex human diseases and traits. Disease is gratefully acknowledged. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


Introduction
Recent data emerging from theoretical models [1,2] and empirical observation through genome-wide association studies (GWAS) (for example [3,4]) demonstrate that hundreds of genetic loci contribute to complex traits in humans. These data prompt two questions: (1) can additional genetic loci be identified by follow-up of the most significantly associated variants after initial GWAS meta-analysis? and (2) can further investigation via genetic fine-mapping refine association signals at established genetic loci? Systematically addressing these two questions should help improve understanding of the genetic architecture of complex traits and their shared genetic determinants, and suggest hypotheses and disease mechanisms that can be tested in functional experiments or model systems [5].
Addressing these two questions requires genotyping thousands of individuals at many genetic markers. For most currently available genotyping technologies, this kind of characterization is costprohibitive. To address this need in the context of type 2 diabetes, coronary artery disease and myocardial infarction, and quantitative traits related to these diseases, we designed the Metabochip, a custom genotyping array that provides accurate and cost-effective genotyping of nearly 200,000 single nucleotide polymorphisms (SNPs) chosen based on GWAS meta-analyses of 23 traits (Table 1). Metabochip SNPs were selected from the catalogs developed by the International HapMap [6] and 1000 Genomes [7] Projects, allowing inclusion of SNPs across a wide range of the allele frequency spectrum. These included 63,450 SNPs to follow-up the top ,5,000 or ,1,000 (see Methods) independent association signals for each of the 23 traits, 122,241 SNPs to fine-map 257 loci which showed genome-wide significant evidence for association with one or more of the 23 traits, and 16,992 SNPs chosen for a variety of other reasons (see Methods and Table 2). In designing the array, we sought to maximize assay success rates as well as the number of variants that could be assayed; Illumina custom arrays include a fixed number of ''beads'' and some sites can be assayed with a single bead while others require two [8].
Here, we describe Metabochip array design, and evaluate performance of the array in common genetic analysis steps, including quality control steps such as genomic control calculations, identification of related individuals, and fine-mapping of known disease susceptibility loci. Our results provide practical guidance to investigators and show that for fine-mapping loci the Metabochip provides much greater resolution than prior GWAS arrays.
We included three design classes of SNPs on the Metabochip ( Table 2): 1. Replication SNPs: ,5,000 (Tier 1) or ,1,000 (Tier 2) SNPs were selected to follow-up the top independent association signals from the largest available GWAS meta-analysis for each of the 23 traits (Supplementary Table S1). 2. Fine-mapping SNPs: SNPs were selected from the catalogs of the International HapMap Project [6] and the August 2009 release of the 1000 Genomes Project [7] to fine-map 257 loci associated at genome-wide significance (P,5610 28 ) in preliminary analyses of one or more of the 23 traits (See Figure 1, Supplementary  Table S2 and S3, and Supplementary Text for details). 3. Other SNPs: These were comprised of independent SNPs for which genome-wide significant associations had been reported for any trait, SNP tags for copy number polymorphisms (CNPs), the MHC region, and the mitochondrial genome, fingerprint SNPs from GWA array products, a set of chromosome X and Y markers for sex verification, and ''wild-card'' SNPs based on consortium-specific hypotheses and interests (for example, based on a known pathway or early deep-sequencing studies). A detailed description of how SNPs were selected in each of these categories can be found in the Supplementary Text [21][22][23][24][25].
In total, 217,695 SNPs were chosen for the array (Table 2). 20,970 SNPs (9.6%) failed during the assay manufacturing process, resulting in 196,725 SNPs available for genotyping. A summary file annotating each Metabochip SNP with ascertainment criteria, SNP assay, a list of unintended duplicate SNPs (Supplementary  Table S4), and reference strand orientation for alleles is provided at http://www.sph.umich.edu/csg/kang/MetaboChip/.

Data Generation and Quality Control (QC)
We evaluated the utility of the Metabochip and accuracy of its genotype calls in three sample sets: (1) 15,896 northern European  individuals from the FUSION, METSIM, HUNT, Tromsø, and Diagen studies [26][27][28][29][30] together with 67 HapMap samples genotyped at least two times each and called using Illumina GenomeStudio software by re-clustering these data; (2) 6,614 Sardinian individuals organized in 1,243 extended families from the SardiNIA study [31,32] called by GenomeStudio software using default cluster data; and (3) 9,715 Nordic individuals from the Malmø Preventive Project, the Scania Diabetes Registry, and the Botnia Study [33][34][35] genotyped using a modified version of the BIRDSEED genotype calling algorithm [36].
We applied standard SNP-and sample-based QC filters based on call rate, Hardy-Weinberg equilibrium deviations, duplicate genotype inconsistencies, and failures of Mendelian inheritance; in the Nordic sample, we also carried out checks based on platespecific characteristics. These filters resulted in final data sets of 163,222 polymorphic SNPs genotyped in 67 HapMap samples, 142,812 polymorphic SNPs genotyped in 6,164 Sardinians, and

Author Summary
Recent genetic studies have identified hundreds of regions of the human genome that contribute to risk for type 2 diabetes, coronary artery disease and myocardial infarction, and to related quantitative traits such as body mass index, glucose and insulin levels, blood lipid levels, and blood pressure. These results motivate two central questions: (1) can further genetic investigation identify additional associated regions?; and (2) can more detailed genetic investigation help us identify the causal variants (or variants more strongly correlated with the causal variants) in the regions identified so far? Addressing these questions requires assaying many genetic variants in DNA samples from thousands of individuals, which is expensive and timeconsuming when done a few SNPs at a time. To facilitate these investigations, we designed the ''Metabochip,'' a custom genotyping array that assays variation in nearly 200,000 sites in the human genome. Here we describe the Metabochip, evaluate its performance in assaying human genetic variation, and describe solutions to methodological challenges commonly encountered in its analysis.

Statistical Analysis Using Metabochip: Genomic Control, PCA, and Kinship Estimation
Since Metabochip SNPs were selected to be associated with our 23 traits of interest, performing genomic control correction [37] requires some care. To select a set of (near)-independent SNPs that are not associated with an analysis trait of interest, we focused on SNPs selected to replicate signals unrelated to the trait of interest (for example, QT interval SNPs for a T2D association analysis), also removing SNPs within 250 kb of SNPs previously associated with the trait of interest, and then LD-pruning the remaining SNPs so that no SNP pair is in strong LD (r 2 ..3).
To estimate kinship coefficients or to correct for population stratification using principal components analysis (PCA) or multidimensional scaling (MDS) covariates, we require SNPs that are not too rare and are not in strong pairwise LD. We found that taking SNPs with MAF..05 and LD-pruning them so that no SNP pair has r 2 ..3 works well for PCA and MDS (data not shown). The same subset of SNPs can be used for pairwise IBD estimation using the maximum-likelihood method of Milligan [38] implemented in PLINK [39] or the variance-components method of Balding and Nichols [40] implemented in EMMAX [41].

Imputation Preparation and Evaluation
We carried out genotype imputation in the Sardinian data. We imputed variants observed in a reference set of 280 Europeans from the August 2010 1000 Genomes Project data into: (a) 6,164 individuals genotyped on the Metabochip [32], (b) 1,097 individuals genotyped on the Affymetrix 6.0 array, and (c) 1,412 individuals genotyped on the Affymetrix 500 K array [42]. We evaluated mean estimated r 2 within fine-mapping regions using minimac ( [43]; www.genome.sph.umich.edu/wiki/minimac), and empirically compared the imputation quality using the published Sanger sequencing data in five fine mapping loci [32]. In addition, we evaluated mean estimated r 2 across different continental populations by leaving one individual out from the 1000 Genomes reference panel and imputing them using markers present in each platform across the fine mapping regions and a 1 Mb window flanking each region. We also compared association power obtained by imputation into GWAS and Metabochip samples in Metabochip fine-mapping regions by comparing LDL cholesterol association evidence in 2,342 of these individuals genotyped using both the Metabochip and one of the Affymetrix arrays.

Evaluation of Array Design and Genotype Quality
Of 217,695 SNPs chosen for the Metabochip across all design categories, 196,725 (90.4%) were successfully manufactured on the array ( Table 2). The 48,846 previously manufactured SNPs had higher success rate (95.4%) than the 168,849 new SNP assays (88.7%). Illumina design score was predictive of the quality of manufactured SNP assays. For example, 25% of SNPs with design score,0.6 failed to produce genotype calls due to poor clustering of the intensity data, compared to 3.1% of SNPs with design score between 0.6 and 1.0 (Supplementary Figure S1).
We evaluated genotype calling accuracy for 67 HapMap samples genotyped multiple times using three different calling strategies: (a) Illumina GenomeStudio with reclustering the intensity data using .15,000 samples; (b) Illumina GenomeStudio based on default clusters provided by Illumina; and (c) GenoSNP [44], which calls genotypes based on a within-sample-betweenmarkers analysis of intensity data rather than a between-samplewithin-marker analysis.
Using GenomeStudio with reclustering, genotype concordance between Metabochip genotypes for duplicate pairs was 99.998% overall and 99.990% for heterozygotes. Comparing Metabochip genotypes to HapMap 3 genotypes for the 59,935 SNPs in common, genotype concordance was 99.93% overall and 99.84% for heterozygotes, similar to the 99.87% Mendelian consistency rate reported in the HapMap3 data [45]. We observed similar concordance rates for these sample sets using the Illumina caller with default clusters (99.93% overall, 99.84% for heterozygotes), or using GenoSNP [44] (99.85% overall, 99.81% for heterozygotes). Genotype concordance for less common variants was slightly lower than for common variants. For example, among the singleton SNPs in the 67 HapMap samples, 98.9% of heterozygous genotypes were concordant with HapMap3 for the two GenomeStudio call sets and 97.8% for the GenoSNP set. Heterozygous genotype concordances for singleton SNPs between duplicate pairs were 99.76%, 99.70%, and 99.83% for the three call sets.

Genotype Imputation within the Metabochip Fine-Mapping Regions
We next investigated accuracy of genotype imputation into the 257 Metabochip fine-mapping regions using the 280 Europeans from 1000 Genomes Project [7] as reference set and the 6,164 individuals in the Sardinian Metabochip sample as target. Figure 3 displays estimated r 2 values in the Metabochip fine-mapping regions as a function of MAF. Also displayed are estimated r 2 values for SNPs in these regions using the 280 European 1000 Genomes project samples as reference set and 1,412 Sardinians genotyped on the Affymetrix 500 K and 1,097 Sardinians genotyped on the Affymetrix 6.0 chips as targets. Imputation accuracy into the Sardinian Metabochip sample is greater in all allele frequency ranges than for the samples genotyped using the GWAS arrays. For example, among SNPs with .02#MAF,.05, mean estimated r 2 for the Affymetrix 500 K, Affymetrix 6.0, and Metabochip samples were .47, .62, and .84, respectively (Figure 4). The improved imputation accuracy for Metabochip compared to GWAS array is primarily due to increased marker density of the Metabochip in these regions. Imputation quality in the Metabochip fine-mapping regions using Metabochip is also improved for non-European individuals compared to imputation using GWAS platforms. Using a leaveone-sample-out approach, we evaluated the average r 2 from the 1000 Genomes reference panel into Affymetrix 500 k, Affymetrix 6.0, and Metabochip. For example, among SNPs with .02,MAF,.05, mean estimated r 2 across European individuals for the chips were .78, .83, and .93, respectively. For individuals with African ancestry, corresponding values were .78, .85, and .94, and for individuals of Asian ancestry, they were .67, .72, and .89 (Supplementary Figure S2). The fact that imputation of rare variants in African ancestry populations is more accurate than in European populations is probably explained by noting that -in the short regions evaluated here -there will be only a limited number of common variant haplotypes in Europeans and, in some cases, these will not effectively tag specific rare variants. In African populations, with a larger variety of rare haplotypes, it is more likely (relative to Europeans) that at least one haplotype will capture rare variants of interest.
In addition, we empirically evaluated the quality of experimentally determined and imputed SNPs within the five fine mapping regions by comparing individual genotypes with those obtained by Sanger sequencing. For 126 SNPs evaluated, the average r 2 in analyses based on the Affymetrix 500 k and 6.0 arrays was .46 and .55, respectively. Analyses based on Metabochip showed average r 2 = .79. Focusing on 48 SNPs that were imputed in all three analyses, the average r 2 was .31 (Affymetrix 500 K), .41 (Affymetrix 6.0), and .57 (Metabochip) (Supplementary Figure S3).

High-Resolution Association Analysis within Metabochip Fine-Mapping Regions
To compare the power and resolution for association testing in the Metabochip fine-mapping regions to that of standard GWAS arrays, we revisited the LDL cholesterol association analysis from the SardiNIA study [32] in 2,342 individuals genotyped for both Metabochip and an Affymetrix (6.0 or 500 k) GWAS chip. Here, we focus on five of the six most strongly associated loci from Willer et al. [46], in and around PCSK9, LDLR, APOE/APOC1/ APOC2, SORT1, and APOB ( Figure 5A-J), all of which were designated for locus fine mapping by the Global Lipids Genetics Consortium.
In the SORT1 and APOB regions, the peak association signals for the two data sets are similar ( Figure 5A-D). For PCSK9, LDLR, and APOE/APOC1/APOC2, Metabochip based analysis resulted in considerably stronger association signals. For PCSK9 and APOE/ APOC1/APOC2, the most strongly associated variants were lowfrequency SNPs (MAF = 1.1% for PCSK9, MAF = 3.4% for APOE) that were directly genotyped on the Metabochip but not on the Affymetrix chips ( Figure 5E-J). Although the signals from common variants are similar, the peak SNPs were not imputed accurately in the Affymetrix data (estimated r 2 = .04 and .08, respectively). Within the LDLR region, there are 165 SNPs in the 1000 Genomes European panel. None of these SNPs are on the Affymetrix chips and only eight could be imputed at estimated r 2 $.3 using the Affymetrix data; the locus is also hard to impute using HapMap 2 as a reference, with the peak association signals corresponding to r 2 of ,. 40. In contrast, 36 of the 165 SNPs were directly genotyped in Metabochip, and 122 were imputed at estimated r 2 $.3. As a result, imputation into the Metabochip data resulted in a substantial association signal (p = 7.3610 26 ), while for the Affymetrix data, p..02 at all markers ( Figure 5I-J). These results demonstrate that dense genotyping may substantially improve imputation accuracy, increasing association power even for common variants.

Performing Standard Statistical Analyses Using Metabochip Genotype Data
We carried out kinship estimation between pairs of individuals and calculated genotype-based principal components for inclusion as covariates in genetic association analysis using all Metabochip SNPs that passed QC, and then using the pruned subset of SNPs described in the Methods section. When using all QC-passing  Table 2) in each of the tabulated minor allele-frequency bins. CNP = copy number polymorphism. doi:10.1371/journal.pgen.1002793.g002 SNPs, estimates of pairwise kinship coefficients in the Sardinia sample had inflated variance (Supplementary Table S5), and kinship coefficient estimates for the Nordic sample calculated using PLINK suggested (incorrectly) that essentially all pairs of individuals were related (Supplementary Figure S4). For each analysis, using the pruned set of SNPs gave sensible results, Because many Metabochip SNPs were included specifically due to prior evidence for association of T2D, CAD/MI and related traits, controlling for potential population stratification in Metabochip analysis requires some care. Not surprisingly, carrying out T2D association analysis in the Nordic sample on all SNPs passing QC without inclusion of genotype-based principal components resulted in a large genomic control inflation factor (l GC = 1.44). Including all SNPs that passed QC to estimate principal components (PCs), and then including those PCs as covariates in the association analysis gave reduced but still substantial inflation (l GC = 1.13). When we instead estimated test statistic inflation based only on the 3,772 LD-pruned QT interval replication SNPs (not expected to associate with T2D) we obtained a genomic control inflation factor near unity (l GC = 1.01).

Assessing Overlap among SNPs across Traits
We were interested whether the replication SNP sets submitted by the GWAS consortia for the different traits showed more or less overlap than expected by chance. To address this question, we counted the number of SNPs in common across pairs of traits, and used simulation to test whether the observed overlaps were different than expected under the null hypothesis of genetic independence of pairs of traits (Supplementary Table S6). Not surprisingly, we observed substantial SNP set overlaps (and greater than expected assuming independence) for multiple pairs of correlated traits, notably SBP and DBP (38% proportion of maximum possible overlap), HDL and TG (17%), and TC and LDL (87%). We also observed substantial genetic overlap (4%) between LDL and SBP, which are nearly uncorrelated traits. Overall, we observed an excess of nominally significant SNP set overlaps, consistent with (but in no way proof of) the hypothesis a shared genetic etiology between these cardiometabolic traits.

Discussion
We designed the Metabochip, a custom genotyping array for replication of the top association signals from the largest available GWAS meta-analysis for 23 T2D and CAD/MI related traits and for fine-mapping 257 genome-wide significant association signals for 15 of these traits ( Table 1). The Metabochip also includes a set of SNPs representing genome-wide significant associations across a range of human traits; SNPs that tag known copy number polymorphisms, the MHC, and mitochondrial variants; X and Y chromosome SNPs for sex verification, fingerprint SNPs for sample tracking, and ''wildcard'' SNPs selected by the participating GWAS consortia ( Table 2). The array has already been genotyped on DNA samples from hundreds of thousands of individuals and preliminary analyses across the contributing GWAS consortia have identified hundreds of new genome-wide association signals (manuscripts being prepared by each of the consortia).
In designing the Metabochip, 90.4% of chosen SNPs were successfully designed and manufactured onto the array, and of these, ,82% passed QC filters in our three example studies, resulting in very complete coverage of variation in our 257 fine-mapping regions. Of course, as time passes and catalogs of SNPs expand, potential shortcomings in coverage should become apparent. Currently, coverage of 1000 Genomes Pilot Study European SNPs in the fine-mapping regions is 82.0% at a tagging threshold of r 2 $.8. Coverage of Phase 1 European SNPs in these regions is 54.5%, and the number increases to 73.7% for SNPs at MAF.0.5%. Using genotype imputation, we can impute 82% of 1000 Genomes Phase 1 European SNPs with MAF.0.5% with estimated r 2 $0.8. The resulting data are of high quality, with 99.99% duplicate consistency in heterozygotes and 99.77% Mendelian consistency in heterozygotes in our studies. Further, Metabochip fine-mapping regions provide an excellent target for genotype imputation from relevant reference sets, and in our experience can provide more complete coverage than provided by standard HapMap-  A key decision in the fine-mapping of any GWAS signal concerns the size of the region where genetic variation will be examined exhaustively. In designing the Metabochip, we focused on relatively small regions surrounding each lead SNP -these included all variants in strong linkage disequilibrium (r 2 ..5) and a small shoulder extending .02 cM beyond that (typically, ,20 kb). This decision was informed by the observation that, in cases where GWAS signals and Mendelian disease loci overlap, they are typically very close together (typically within ,10 kb of each other and nearly always within ,100 kb; see [4] for a discussion of the issue), although there are exceptions to this rule (see [47], for example).
Within each fine-mapping region, we selected variants identified by the HapMap consortium and early analyses of the 1000 Genomes Consortium data. The 1000 Genomes Project and other sequence based catalogs of genetic variation are now more extensive that at the time of array design, but (as noted above) our analyses show that the SNPs selected for inclusion in the Metabochip form a useful reagent for genotyping imputationnot only for the imputation of newly discovered SNPs in the finemapping regions (see above) but also for the imputation of other types of variants, such as indel polymorphisms, that have become part of newer 1000 Genomes Project analyses (unpublished data).
Several other design choices for Metabochip were to some degree arbitrary: which traits to include; balance in numbers of SNPs for replication, fine mapping, and other purposes; and how to prioritize among SNPs available for each purpose. Were we to design a similar chip now, we would take advantage of the now available more extensive and deeply annotated SNP catalogs. In addition, we would likely include a set of randomly ascertained SNPs to facilitate analysis that control for population structure and other artifacts. Finally, with empirical evidence from this and other projects on the relationship between SNP design score and empirical probability of successful design, we would likely replace design score by probability of successful design. This approach would likely result in even higher call rates.
Because Metabochip SNPs are highly enriched for traitassociated SNPs and .60% are clustered in the ,1.5% of the genome that comprises the fine-mapping regions, Metabochip genotype data present some challenges to standard analyses such as relationship estimation, principal components analysis, and genomic control determination. However, as we demonstrated, these challenges can be overcome by focusing on replication SNPs expected to be unrelated to the trait of interest. An alternative approach is to use SNPs that were not associated with the trait(s) of interest in the corresponding GWAS (for example, p-value..50 for all such traits) and then to LD-prune the resulting set of SNPs to identify a near-independent set. An alternative that is also worthy of investigation in the analysis of case-control samples is the application of principal component factor loadings derived from a controls-only analysis to the combined sample of cases and controls. When this last alternative is considered, it is important to check that PCA axes derived from controls represent all relevant ancestries present in cases. The design of the array, focused on replication and fine-mapping and selecting SNPs from early releases of the HapMap and 1000 Genomes Projects, resulted in a highly non-random ascertainment of SNPs. Thus, we cannot recommend use of Metabochip SNPs for population genetic analyses that rely on unbiased, and/or comprehensive ascertainment schemes for SNPs.
The need for follow-up genotyping is a frequent requirement of GWAS and sequencing studies of complex human traits. Approaching array design in a coordinated fashion across related studies and traits can be particularly cost-effective, since per array costs often drop dramatically with increasing numbers of individuals to be genotyped, and (given sufficient numbers of individuals) may increase only modestly with increasing numbers of SNPs. For example, a custom chip designed to genotype the ,22,000 DIAGRAM-selected type 2 diabetes Metabochip SNPs in the ,80,000 individuals genotyped on Metabochip by the DIAGRAM consortium studies would have cost ,$55 compared to the Metabochip cost of $39, delivering only 1/9 as many genotypes at .40% greater cost. Furthermore, examining the association between SNPs tentatively associated with one trait for other related traits can also be informative, highlighting pleiotropy across related traits and helping discover new association signals; for example, two of the ten novel type 2 diabetes loci identified to date by Metabochip analysis by the DIAGRAM consortium were placed on Metabochip for other traits [48]. In the case of the Metabochip, which is less expensive than many smaller trait specific arrays, this opportunity to collect more information and investigate the effects of SNPs associated with other traits actually comes with reduced costs (compared to trait specific arrays), although with the need to organize across multiple consortia and to share the number of SNPs that can be cost-effectively genotyped. The ''Immunochip'' [49] follows this same paradigm and supports genotyping of ,200,000 SNPs identified on the basis of GWAS meta-analyses for immunological disorders, while the recently designed ''exome chip'' (Benjamin Neale, Gonçalo Abecasis, personal communication) supports genotyping of ,250,000 exonic SNPs identified via large-scale exome sequencing studies totaling .12,000 individuals. These and other similar array products represent valuable tools in ongoing efforts to understand the genetic architecture of complex human traits.    Text S1 Technical details of SNP selection criteria. (DOC)