As most of the heritability of complex traits is attributed to common and low frequency genetic variants, imputing them by combining genotyping chips and large sequenced reference panels is the most cost-effective approach to discover the genetic basis of these traits. Association summary statistics from genome-wide meta-analyses are available for hundreds of traits. Updating these to ever-increasing reference panels is very cumbersome as it requires reimputation of the genetic data, rerunning the association scan, and meta-analysing the results. A much more efficient method is to directly impute the summary statistics, termed as summary statistics imputation, which we improved to accommodate variable sample size across SNVs. Its performance relative to genotype imputation and practical utility has not yet been fully investigated. To this end, we compared the two approaches on real (genotyped and imputed) data from 120K samples from the UK Biobank and show that, genotype imputation boasts a 3- to 5-fold lower root-mean-square error, and better distinguishes true associations from null ones: We observed the largest differences in power for variants with low minor allele frequency and low imputation quality. For fixed false positive rates of 0. 001,0. 01,0. 05, using summary statistics imputation yielded a decrease in statistical power by 9,43 and 35%, respectively. To test its capacity to discover novel associations, we applied summary statistics imputation to the GIANT height meta-analysis summary statistics covering HapMap variants, and identified 34 novel loci, 19 of which replicated using data in the UK Biobank. Additionally, we successfully replicated 55 out of the 111 variants published in an exome chip study. Our study demonstrates that summary statistics imputation is a very efficient and cost-effective way to identify and fine-map trait-associated loci. Moreover, the ability to impute summary statistics is important for follow-up analyses, such as Mendelian randomisation or LD-score regression.
Genome-wide association studies (GWASs) have been successfully applied to reveal genetic markers associated with hundreds of traits and diseases. The genotyping arrays used in these studies only interrogate a small proportion of the genome and are therefore typically unable to pinpoint the causal variant. Such arrays have been designed to be cost-effective and include only a set of tag single nucleotide variants (SNVs) that allow the inference of many other unmeasured markers. To date, thousands of individuals have been sequenced [1,2] to provide high resolution haplotypes for genotype imputation tools such as IMPUTE and minimac [3,4], which are able to infer sequence variants with ever-increasing accuracy as the reference haplotype set grows. Downstream analyses such as Mendelian randomisation [5], approximate conditional analysis [6], heritability estimation [7], and enrichment analysis using high resolution annotation (such as DHS) [8] often require genome-wide association results at the highest possible genomic resolution. Summary statistics imputation [9] has been proposed as a solution that only requires summary statistics and the linkage disequilibrium (LD) information estimated from the latest sequencing panel to directly impute up-to-date meta-analysis summary statistics [10]. Because summary statistics imputation uses summarised data as input, it is not bounded to privacy restrictions related to the use of individual data. Another advantage is its substantially lower computation time compared to genotype imputation. For example, for imputation of the UK Biobank data, it is about 500 times faster (4200 vs 8. 3 CPU days comparing Minimac [4] to our SSIMP software [11]). This study compares summary statistics imputation directly to genotype imputation and focuses on its practical advantages using real data. In particular, we evaluated two experiments: 1) we ran a GWAS on both simulated traits and human height using data from 120′086 individuals from the UK Biobank and compared the performances of summary statistics imputation and genotype imputation, using direct genotyping/sequencing as gold standard; 2) we imputed association summary statistics from a HapMap-based GWAS study [12] using the UK10K reference panel to explore new potential height-associated variants which we validated using results from Marouli et al. [13] and the UK Biobank height GWAS (n = 336′474). We extended summary statistics imputation [9,14] which yields increased imputation accuracy by accounting for variable sample sizes. For all applications presented in this manuscript we are using this improved version of summary statistics imputation.
By combining summary statistics for a set of variants and the fine-scale LD structure in the same region, we can estimate summary statistics of new, untyped variants at the same locus. We assume a set of univariate effect size estimates ai are available for SNVs i = 1, …, I from a linear regression between a continuous phenotype y and the corresponding genotype gi measured in N individuals. Without loss of generality we assume that both vectors are normalised to have zero mean and unit variance. Thus a i = (g i) ′ · y N and a= (a1, a2, …, aI) ′∼N (α, Σ). Σ represents the pairwise covariance matrix of effect sizes of all i = 1, …, I SNVs. To estimate the univariate effect size αu of an untyped SNV u in the same sample, one can use the conditional expectation of a multivariate normal distribution. The conditional mean of the effect of SNV u can be expressed using the effect size estimates of the tag SNVs [9,15]: a ^ u = a u | M = α u + Σ u M Σ M M - 1 (a - α), (1) where M is a vector of so-called tag SNVs, Σ u M represents the covariance between SNV u and all M markers and Σ M M represents the covariance between all M markers. We assume that estimates for the two covariances are available from an external reference panel with n individuals and denote them s = Σ ^ M u, S = Σ ^ M M. The corresponding correlation matrices are γ and Γ, with c = N ⋅ s and C = N ⋅ S being the estimates for the correlation matrices. Further, by assuming that SNV u and the trait are independent conditioned on the M markers, i. e. α u - Σ u M Σ M M - 1 α = 0, Eq (1) becomes a ^ u = a u | M = s ′ S - 1 a = c ′ C - 1 a (2) One can also choose to impute the Z-statistic instead, as derived by Pasaniuc et al. [9]: z ^ u | M = c ′ C - 1 z (3) with z = a N, when the effect size is small (as is the case in typical GWAS). Similar to Pasaniuc et al. [9], we chose M to include all measured variants within at least 250 Kb of SNV u. To speed up the computation when imputing SNVs genome-wide, we apply a windowing strategy, where SNVs within a 1 Mb window are imputed simultaneously using the same set of M tag SNVs the 1 Mb window plus 250 Kb flanking regions on each side. We simulated genetic data on 25’000 individuals was used. In brief, we used data from the five European subpopulations CEU, GBR, FIN, TSI and IBR of the 1000 Genomes reference panel [1]. We chose to up-sample chromosome 15 using HAPGEN2 [29] to 5′000 individuals for each subpopulation, yielding a total of 25′000 individuals. Of these, half of the data was used to estimate the LD structure C and the other half to simulate the association study with an in silico phenotype. The simulation procedure is described in more detail in S1 Appendix. Forty regions were selected with one non-HapMap causal variant in each and all HapMap SNVs were used as tag SNVs. Sample size distributions were drawn from two published GWAS studies (on HDL [30] and T2D [31]). Missingness was assigned at random positions while respecting the missingness correlation parameter θmiss, with zero value reflecting missingness at random and one corresponding to the maximum possible sample overlap between SNVs. To estimate LD structure in C and c (Eq (2) ) we used 3′781 individuals from UK10K data [32,33], a reference panel of British ancestry that combines the TWINSUK and ALSPAC cohorts. All analysis was performed with R-3. 2. 5 [34] programming language, except GWAS summary statistics computation for UK Biobank genotype and genotype imputed data, for which SNPTEST-5. 2 [35] was used. For summary statistics imputation we used SSIMP [11].
The conventional estimate of the standardised effect of a SNV u, a ^ u (c o n v), (Eq (2) ) is unbiased, under certain assumptions, but can have large variance when there is variation in the sample sizes recorded in N M. In this section, we used upsampled 1000 Genomes data [1] and simulated phenotype with known standardised effect α and various missingness design. We compare the MSE of the conventional estimation to the MSE of two other estimators, Eq (13) using D (dep) and D (ind), derived in the method section. In general, the size of the overlap is unknown and we recommend using the assumption of maximum dependence (D (dep) ) as it is the most conservative assumption. An alternative is to assume randomly distributed missingness (D (ind) ). Most pairs of SNVs in GIANT attain close to the maximum possible missingness-overlap (S10 Fig) and therefore this assumption is not overly-conservative. The results in Fig 3 demonstrate that the conventional method has the largest MSE across all the simulation parameters tested. Where the variance in sample size is very large (top row of Fig 3), the true correlation is often very close to zero. Both of our methods effectively make this same (correct) assumption of low correlation and therefore they both perform equally well. Where the variation in sample size is less extreme, as in the simulations on the bottom row of Fig 3, there is less shrinkage of correlation and the simulated missingness correlation becomes more relevant. Where the simulated data has the maximum possible missingness correlation (on the right hand side of the subplots in Fig 3), i. e. the sample overlap between each pair of SNVs is as large as possible given their two sample sizes, D (dep) performs better (as expected). With lower overlap (first column) D (ind) performs better. By having two types of genetic data at hand, genotype and imputed genotype data, we were able to compare summary statistics of 37′467 typed SNVs resulting from (1) associations calculated from original genotype data (ground truth); (2) associations calculated from imputed genotype data (genotype imputation) and (3) associations imputed from summary statistics calculated using genotype data (Fig 1). For our analysis, we defined 706 genomic regions in total, among which 535 contain SNVs associated with height [12,13], while the remaining 171 regions were selected to be free of any known height associated SNVs. We examined imputation results for different SNV categories. These were grouped based on (i) their association status (being correlated with the causal SNV vs. null SNVs) with the lead SNV of each of the 535 height-associated regions (6′080 variants were correlated, 31′567 were not); (ii) frequency (MAF: 1% < low-frequency ≤ 5% < common; 13′857 and 23′790 variants, respectively); and (iii) imputation quality based on summary statistics imputation (r ^ pred, adj 2: low ≤ 0. 3 < medium ≤ 0. 7 < high; 724,9′792, and 27′131 variants, respectively). S1 and S2 Figs show the distribution of SNV counts in each of these twelve subgroups. We term the 6′080 SNVs correlated with a height-associated lead SNV as associated SNVs. Conversely, we refer to the 31′567 SNVs that are not correlated with any height-associated lead SNV as null SNVs. For both, null and associated SNV groups, the largest group of analysed variants were common and well-imputed (S1 Fig). The fraction of SNVs with low quality imputation increases with lower minor allele frequency (S2 Fig). However, the number of rare variants (MAF < 1%) were too small (2′411 variants, among these only 13 associated variants), similar to the number of badly-imputed SNVs (724 variants, among these only one associated variant) to draw meaningful conclusions and hence we limited our analysis to common and low-frequency, and medium- and well-imputed variants. We focused on two aspects of the imputation results. First, we compared how summary statistics imputation and genotype imputation perform relative to the ground truth (direct genotyping). For this we used four measures: the root mean squared error (RMSE), bias, the linear regression slope, and the correlation. Second, we calculated power and false positive rate for genotype imputation and summary statistics imputation directly. While previous studies have examined the role of (common) HapMap variants for height [12,36], the impact of rare coding variants could not be investigated until bespoke genotyping chips (interrogating low-frequency and rare coding variants) were designed to address this question in a cost-effective manner. Such an exome chip based study was conducted by the GIANT consortium in 381′000 individuals and revealed 120 height-associated loci, of which 83 loci were rare or low-frequency [13]. These association results enabled us to compare the usefulness of imputation-based inference with direct genotyping done in Wood et al. [12], since the two studies are highly comparable in terms of ancestry composition and statistical analysis, evidenced by S6 Fig confirming very high concordance between summary statistics for the subset of 2′601 SNVs correlated to a height-associated variant which were available in both studies.
Imputation accuracy is affected by the varying sample size across tag SNVs. If two SNVs were observed in two different samples, the correlation between the summary statistics will decrease with the number of individuals in common between the two samples. Our approach addresses this problem by shrinking the correlation matrix according to sample size overlap. We present two ways of estimating this overlap: D (ind) for independent missingness, which is randomly distributed; and D (dep) for dependent missingness, which is highly correlated. To evaluate the performance ot these two methods we simulated data with two different distributions of missingness (narrow or wide range of sample sizes) and varying correlation in missingness between variants (from completely random to maximal overlap, Fig 3). We then compared the performances of conventional summary statistics imputation and our proposed dependent (D (dep) ) and independent (D (ind) ) approaches. Overall, replacing C and c with D and d yields a lower RMSE. Furthermore, we note that the dependent approach has lower RMSE when the sample size variance is low and the missingness correlation approaches one. S15 Fig shows the comparison between the conventional estimation and using D (dep) for imputing GIANT height association summary statistics. Ideally, for any pair of SNVs that are in LD with each other, we would know the exact number of individuals that are in the overlap, i. e. the number of individuals for which both SNVs were genotyped. Using the individual study missingness and sample sizes from the Genetic Investigation of ANthropometric Traits (GIANT) consortium, we demonstrate in Fig. S10 Fig that the size of the overlap is generally larger than would be the case under a strict ‘missing independently at random’ assumption. Furthermore, the correlation of missingness is typically positive (N k ∩ l > N k N l N m a x) and often approaches the maximum possible overlap (Nk∩l = min (Nk, Nl) ). The reason for this is that SNPs are either entirely missing from a study or being available for all study participants depending on its genotyping chip or imputation panel, which induces positive missingness correlation between markers. We compared summary statistics imputation and genotype imputation by using individual-level data from the UK Biobank. In general, imputation using summary statistics imputation leads to a larger RMSE than genotype imputation in all twelve SNV subgroups investigated (Fig 6). Among associated SNVs, summary statistics imputation performs similar to genotype imputation for well-imputed SNVs, but shows a trend for underestimation of the Z-statistics and lower correlation with the true effect size for medium-imputed SNVs (Fig 4). Conversely, genotype imputation has more consistent results for most of the twelve SNV subgroups (Figs 4 and 5), that is reflected in a correlation close to one between Z-statistics from genotype data and genotype imputation data. When investigating power and FPR for both methods (Fig 7) we observe that for a given significance threshold, summary statistics imputation has lower power compared to genotype imputation, an effect that is amplified for SNVs with lower imputation quality (r ^ pred, adj 2 ≤ 0. 7) and lower MAF (MAF ≤ 5%). Ultimately, the underestimation of imputed Z-statistics with summary statistics imputation leads to a lower type I error. This effect is amplified for SNV groups with lower imputation quality (r ^ pred, adj 2 < 1). For associated SNVs with r ^ pred, adj 2 < 1 we expect an underestimation for associated SNVs due to the fact that we are imputing summary statistics under the null model, whereas for null SNVs with r ^ pred, adj 2 < 1 we expect an underestimation due to decreased variance of the summary statistics imputation estimation. Ideally, for an unbiased estimation of causal and null SNVs, the imputed Z-statistics (Eq (2) ) should be divided by r ^ 2. However, as the imputation quality r ^ pred, adj 2 is noisily estimated from small reference panels (discussed below) and it is not guaranteed that the SNV we impute is causal, we risk to overestimate the summary statistics of associated SNVs. This is the reason why refrain from doing so. S9 Fig shows the P-value distribution of summary statistics imputation for null SNVs with an accumulation of low P-values for well-imputed SNVs and an accumulation of high P-values for badly-imputed SNVs. We think that two factors are in play here. First, mostly due to polygenicity, the genomic lambda for height is λGC = 1. 94, therefore we expect even seemingly null variants to show inflation. Second, for null SNVs, the sample variance of the imputed Z-statistics should be proportional to the average imputation quality. We calculated for each of the null SNV subgroups the ratio between the sample variance for Z-statistics from summary statistics imputation and the sample variance for Z-statistics from genotype data. For common null SNVs we observe a ratio that gradually decreases with imputation quality (0. 86 for perfectly-, 0. 61 for medium- and 0. 32 for badly imputed SNVs). For low-frequency null variants the ratio is up to 0. 6 lower (0. 80 for perfectly-, 0. 54 for medium- and 0. 30 for badly imputed SNVs). The inflation for well-imputed SNVs can be explained by the genomic lambda, while for badly-imputed SNVs it is aggravated by the underestimated standard error. Because the number of associated SNVs with MAF < 1% was too low (13 variants) to draw any meaningful conclusions, we refrained from analysing this MAF group. One other reason to exclude rare variants from this analysis is, that the reference panel used (UK10K) contains 3′871 individuals and therefore estimations for LD of rare variants are unreliable and rare variants can (in theory) only be covered down to MAF = 1/ (2 ⋅ 3′871). We believe improving summary statistics imputation for rare variants will require not only larger reference panels to allow estimation of LD of rare variants, but also methods which would allow non-linear tagging of variants. It should be kept in mind that, just like for genotype imputation, even with very large reference panels, one will not be able to impute variants with extremely rare allele counts. To investigate these SNVs full genome sequencing is indispensable [42]. We find that our imputation quality measure r ^ pred, adj 2 is conservative and probably underestimates the true imputation quality (S4 Fig). To calculate the imputation quality r ^ pred, adj 2, we need—similar to imputing summary statistics in Eq (2) —to compute correlation matrices c and C estimated from a reference panel (Eq (8) ) and therefore encounter similar challenges as summary statistic imputation itself due to difficulties of reliable LD estimation. The discrepancy in imputation quality metric between summary statistics imputation and genotype imputation (S4 Fig) can be explained by the fact that: (1) genotyped variants that were imputed too, were also used for phasing, (2) it is indeed more difficult to impute summary statistics using summary statistics imputation, and therefore the imputation quality is shifted towards zero, and (3) r ^ pred, adj 2 is an estimation that can either be erroneous due to choosing the wrong reference panel (and therefore r ^ pred, adj 2 does not represent the true imputation quality) or it can be imprecise due to small sample size of the reference panel. For example, UK10K contains 3′871 individuals and is too small to precisely estimate these matrices (the standard error for a correlation estimated from n = 3′871 is 0. 016), which becomes problematic in cases of low correlation. As a showcase of the utility of summary statistics imputation we imputed Wood et al. [12] to higher genomic resolution (limited to variants with MAF ≥ 0. 1% as well as 111 previously reported exome variants) [13], then selected imputed variants that act independently from all variants reported in Wood et al. and from each HapMap SNP, we then replicated these using (independent) UK Biobank data. While Wood et al. [12] is the largest height study to date in terms of number of markers (covering HapMap variants in 253′288 individuals), Marouli et al. [13] exceeds their sample size by more than 100′000 individuals, but is limited to 241′419 exome variants. The similarity between both GIANT studies made the exome chip study ideal for replication. We chose the UK Biobank as a second replication dataset, despite its limitation to individuals of British ancestry, as it covers more variants than the exome chip study. We identify 35 regions, of which one had already been identified in the recent GIANT height exome chip study (rs28929474) and 19 replicated in UK Biobank (at α = 0. 05/35 level). Two candidate loci (#2, #3 in Table 2) that replicate in UK Biobank have borderline significant HapMap signals in close proximity (P-value between 10−6 and 10−8 in [12]) and were therefore not reported in the study in 2014. The 15 non-replicating candidate loci were on average on a lower allele frequency spectrum (ten are rare, three are low-frequency variants, and two are common). Allele frequency was higher among the 20 replicating candidate variants (19 were common and one a low-frequency variant). We also ran an additional approximate conditional analysis, where we conditioned each of the 35 variants onto their neighbouring HapMap SNP (one-by-one). The resulting maximum conditional P-value per locus, is provided as an additional column S1 Table. Correcting for the testing of 529 windows (α = 0. 05/529) we find evidence that 18 of the 35 variants are not only independent from all [12] reported SNPs, but also of each HapMap variant too. For replication of summary statistics from European individuals we use the UK Biobank, which represents only a subset of all European ancestries and is genotype-imputed (instead of sequenced), but on the other hand provides a reliable resource due to its sample size. Furthermore, in UK Biobank, genotype imputation done for genotyped variants can only partially be compared to genotype imputation for untyped variants, as genotyped variants were used for phasing (therefore genotype imputation of genotyped variants is easier and leads imputation qualities close to one, S4 Fig). Due to the small number of height-associated rare variants (13) we can not draw meaningful conclusions for this group and hence avoided their analysis. The choice of the reference panel to conduct summary statistics imputation depends on the fine balance between maximising the sample size of the reference panel (which determines the error in estimated LD) and matching the population diversity of the conducted GWAS. At the first glance, 1000 Genomes reference panel could have been used to appropriately match GIANT allele frequencies, however, the 8-fold higher sample size of UK10K panel offers a larger benefit, ultimately reducing the RMSE [43]. For the simulation study comparing standard summary statistics imputation to our method taking into account variable missingness, we used an upsampling technique called HAPGEN2 [29], which limits the lower bound of the global allele frequency to 1/ (2 ⋅ 503). Furthermore, the outcome used for the simulated GWAS is based on one causal variant with an explained variance of 0. 02, therefore it might not be fully representative for a polygenic phenotype with more than one causal variant. The summary statistics imputation method itself has several limitations too. Due to the size of publicly available sequenced reference panels we can not explore the performance of rare variants (MAF < 1%). The imputation of summary statistics of an untyped SNV is essentially the linear combination of the summary statistics of the tag SNVs (Eq (2) ). Such a model cannot capture non-linear dependence between tag- and target SNVs [10], which is often the case for rare variants [44,45]. In contrast, genotype imputation is able to capture such non-linear relationships by estimating the underlying haplotypes (a non-linear combination of tagging alleles). Furthermore, in case of genotype imputation it is sufficient that the relevant haplotypes are present in the reference panel, but the overall allele frequency does not need to match the GWAS allele frequency. Summary statistics imputation relies on fine tuning of parameters, such as shrinkage of the correlation matrix. Any λ > 0 will make the correlation matrix invertible, but a stronger shrinkage can compensate for estimation error. We hypothesised that optimal shrinkage depends on local LD structure, and sought to optimise λ for each genomic region using the effect sizes of tag SNVs as training data set in a leave-one-out fashion. When looking at null variants, however, maximum shrinkage (λ = 1) usually leads to the smallest RMSE. Therefore, when looking at a region with a mixture of null and associated SNVs, the selected λ will be shifted towards 1 and shrink the estimation of associated SNVs towards 0, which is not ideal. The imputation quality metric r ^ pred, adj 2 tends to be inaccurate in case of small reference panels. The metric is commonly estimated as the total explained variance of a linear model given the reference panel, where the unmeasured SNV is regressed onto all measured markers in the reference panel (Eq (7) ). We noticed that for reference panel sizes smaller than 1000 individuals, the conventional estimation of imputation quality in Eq (7) is biased towards overestimation. We extend the existing imputation quality (Eq (7) ) by accounting for sample size and the effective number of variants (Eq (8) ). The most accurate imputation quality estimations are obtained using an out-of-sample prediction after model selection by fitting a ridge regression model for each unmeasured SNV (r 2 ^ ridge). However, due to the computational complexity, the calculation takes longer than the actual imputation. We provide a more detailed analysis in S2 Appendix.
Genome-wide association studies (GWASs) quantify the effect of genetic variants and traits, such as height. Such estimates are called association summary statistics and are typically publicly shared through publication. Typically, GWASs are carried out by genotyping ∼ 500′000 SNVs for each individual which are then combined with sequenced reference panels to infer untyped SNVs in each’ individuals genome. This process of genotype imputation is resource intensive and can therefore be a limitation when combining many GWASs. An alternative approach is to bypass the use of individual data and directly impute summary statistics. In our work we compare the performance of summary statistics imputation to genotype imputation. We observe that genotype imputation shows a 3- to 5-fold lower RMSE compared to summary statistics imputation, as well as a better capability to distinguish true associations from null results. Furthermore, we demonstrate the potential of summary statistics imputation by presenting 34 novel height-associated loci, 19 of which were confirmed in UK Biobank. Our study demonstrates that given current reference panels, summary statistics imputation is a very efficient and cost-effective way to identify common or low-frequency trait-associated loci.
genome-wide association studies
variant genotypes
sociology
social sciences
alleles
genetic mapping
mathematics
statistics (mathematics)
genome analysis
molecular biology techniques
genotyping
research and analysis methods
mathematical and statistical techniques
statistical methods
consortia
molecular biology
genetic loci
heredity
meta-analysis
genetics
biology and life sciences
physical sciences
genomics
genomics statistics
computational biology
human genetics
Evaluation and application of summary statistic imputation to discover new height-associated loci