Improved imputation of low-frequency and rare variants using the UK10K haplotype reference panel.
Free full text in Europe PMC
Abstract
Free full text
Improved imputation of low-frequency and rare variants using the UK10K haplotype reference panel
Abstract
Imputing genotypes from reference panels created by whole-genome sequencing (WGS) provides a cost-effective strategy for augmenting the single-nucleotide polymorphism (SNP) content of genome-wide arrays. The UK10K Cohorts project has generated a data set of 3,781 whole genomes sequenced at low depth (average 7x), aiming to exhaustively characterize genetic variation down to 0.1% minor allele frequency in the British population. Here we demonstrate the value of this resource for improving imputation accuracy at rare and low-frequency variants in both a UK and an Italian population. We show that large increases in imputation accuracy can be achieved by re-phasing WGS reference panels after initial genotype calling. We also present a method for combining WGS panels to improve variant coverage and downstream imputation accuracy, which we illustrate by integrating 7,562 WGS haplotypes from the UK10K project with 2,184 haplotypes from the 1000 Genomes Project. Finally, we introduce a novel approximation that maintains speed without sacrificing imputation accuracy for rare variants.
Statistical inference of missing genotypes (imputation), where genotyped markers from SNP arrays are used to impute unobserved genotypes from haplotype panels such as the HapMap data, has been instrumental to the discovery of thousands of complex trait loci in meta-analyses of genome-wide association studies (GWAS)1,2. Whole-genome sequencing (WGS) provides near-complete characterization of genetic variation, but it is still prohibitive for researchers to conduct WGS on the large number of samples that are needed to study phenotypic associations of low-frequency and rare genetic variants (minor allele frequency (MAF) <1–5% and <1% respectively). Recently, the 1000 Genomes Project (1000GP) has provided phased haplotypes for more than a thousand samples from diverse worldwide populations, thereby boosting variant coverage and imputation quality, particularly for variants with MAFs of 1–5% (ref. 3). Imputation using this large reference panel has been made computationally efficient by pre-phasing of GWAS samples4 and approximations that select a subset of reference haplotypes5.
Here we describe a novel WGS imputation panel comprising 3,781 samples from the UK10K Cohorts project6. We show that this reference panel greatly increases accuracy and coverage of low-frequency variants relative to a panel of 1,092 individuals from the 1000GP. In addition, we show that imputation accuracy can improve substantially when reference haplotypes are re-phased after initial WGS genotype calling. We present a practical solution for combining imputation reference panels to increase variant coverage, and we introduce a new approximation that maintains the speed of existing approximations while achieving higher accuracy.
Results
The UK10K imputation panel
The UK10K Cohorts Project6 includes two population samples from the UK (http://www.uk10k.org/studies/cohorts.html). The TwinsUK registry comprises unselected, mostly female volunteers ascertained from the general population through national media campaigns in the UK7. The Avon Longitudinal Study of Parents and Children (ALSPAC) is a population-based birth cohort study that recruited >13,000 pregnant women resident in Bristol (formerly Avon), UK8. A total of 1,990 individuals from TwinsUK and 2,040 individuals from ALSPAC were consented for sequencing. Variant sites and genotype likelihoods were called using SAMtools9, and genotypes were refined and phased using Beagle10, following similar procedures to the 1000GP (Methods)3. After QC, 45,492,035 variant sites (42,001,210 single-nucleotide variants and 3,490,825 insertion/deletions (INDELs)) were retained (Table 1) in 1,854 and 1,927 individuals in the TwinsUK and ALSPAC panels, respectively. We downloaded phased haplotypes from 1000GP (Phase 1 integrated v3), which include a total of 39,527,072 sites. We developed new software functionality for merging haplotype reference panels (Supplementary Note 1 and Supplementary Fig. 1). For imputation using the merged panel, here we removed multi-allelic sites and further excluded variants seen only in 1000GP or seen only once in the combined 1000GP+UK10K data set (singletons, see footnote of Table 1 for details). The choice of removing 1000GP-only and singleton sites was designed to specifically evaluate the impact of the increased European-ancestry panel in UK10K vis-à-vis the smaller 1000GP EUR panel. A total of 26,032,603 sites were retained for the imputation reference panel of UK10K panel, and 32,449,428 sites for the imputation reference panel of 1000GP. Given that 16,122,337 sites exist in both panels, combining the two reference panels results in a total of 42,359,694 sites. Overall, 5,775,752 (35.8%) of the overlapping sites had frequencies >5% and another 2,451,738 (15.2%) had frequencies between 1 and 5% in the UK10K sample.
Table 1
| UK10K | 1000GP(Phase 1 v3) | Combined | Overlap | |
|---|---|---|---|---|
| N samples (% European) | 3,781 (100%) | 1,092 (34.7%) | 4,873 | — |
| N total sites in final release | 45,492,035 | 39,527,072 | — | |
| N total sites after filtering* | 26,032,603 | 32,449,428 | 42,359,694 | 16,122,337 |
| Autosome SNPs | 23,411,635 | 29,797,220 | 38,238,102 | 14,970,753 |
| Autosome INDELs | 1,698,262 | 1,370,819 | 2,407,858 | 661,223 |
| Chr X SNPs | 858,380 | 1,223,328 | 1,612,230 | 469,478 |
| Chr X INDELs | 64,326 | 58,061 | 101,504 | 20,883 |
*For UK10K, the following sites were excluded: 18,180,633 singletons that do not exist in 1000GP, 1,064,168 multi-allelic sites and 214,631 mis-matched alleles sites. For 1000GP, the following sites were excluded: 7,053,246 singletons that do not exist in UK10K, 23,932 sites with a SNP and an INDEL at the same position and 443 within large structural deletions. The bold indicates that these four categories of variants are subsets of the N total sites after filtering.
Imputation evaluation of UK10K versus 1000GP
As a first assessment of the UK10K reference panel, we performed a leave-one-out cross-validation on a pseudo-GWAS of UK ancestry, corresponding to a sub-sample of 1,000 individuals from the UK10K WGS data set (500 from TwinsUK and 500 from ALSPAC). For this experiment, we removed each sample from the reference panel in turn, selected 13,413 sites on chromosome 20 from the Illumina 610
k bead chip (pseudo-GWAS panel), and imputed all other sites on this chromosome from a given reference panel. We conducted the imputation with three haplotype reference panels: the 1000GP Phase 1 panel, the ‘original' UK10K panel produced by initial genotype refinement and haplotyping with BEAGLE, and a ‘re-phased' UK10K panel that was generated by using SHAPEIT v2 (ref. 11) to estimate haplotypes from the BEAGLE genotypes (Supplementary Fig. 2). The accuracy of imputed variants was calculated as the squared Pearson correlation coefficient (r2) between imputed genotype dosages in (0–2) and masked sequence genotypes in (0,1,2). The results were stratified into non-overlapping MAF bins for plotting.
The results of this experiment are shown in Fig. 1a, which focuses on variants with MAF<5%. The corresponding plot for all MAF is shown in Supplementary Fig. 3. Both UK10K reference panels (blue dotted and solid lines) produced higher accuracy than the 1000GP panel (black line), with greater gains at lower frequencies. These trends were expected due to the larger sample size and better ancestry matching of the UK10K reference panel to the pseudo-GWAS data. Notably, the UK10K reference panel yielded much higher imputation accuracy after re-phasing with SHAPEIT v2 (solid versus dotted blue lines): the mean r2 at low frequencies increased by >0.1 (20%) after re-phasing, which implies a substantial boost in the power to detect associations. A large imputation panel is a resource that can inform a variety of association studies, so these results suggest that taking the time to improve a WGS panel's haplotype quality could have substantial downstream benefits.
(a) Imputation accuracy in the UK10K pseudo-GWAS test panel using reference panels from 1000GP (black) and UK10K (blue). The ‘original' UK10K reference panel (dotted blue line) was produced by standard genotype refinement of low-coverage sequencing data, whereas the ‘rephased' reference panel (solid blue line) was produced by running SHAPEIT v2 on the genotypes called by BEAGLE to improve haplotype accuracy. (b) Number of imputed variants in UK10K pseudo-GWAS panel as a function of predicted minor allele frequency in the study cohort (x-axis), expected imputation r2 (density of shading), and reference panel: 1000GP (black), UK10K (blue), or combined UK10K and 1000GP (red). Confidently imputed variants are shown in the bottom segment of each bar for easy comparison. Note that expected r2 tends to be larger than true r2. (c) As in b, but using the INCIPE cohort (representative of the general Italian population) as a pseudo-GWAS panel. (d) Imputation accuracy in the INCIPE pseudo-GWAS panel using the UK10K reference panel and different imputation approximations. Results are provided for a run that used all reference haplotypes with no approximation (blue solid line), a run that used an established Hamming distance approximation (orange solid line), and a run that used a new tract sharing approximation (orange dashed line).
Evaluation of combining two reference panels
It is becoming increasingly common for investigators to conduct their own WGS of particular study populations, and a natural goal is to combine these data sets with publicly available reference panels (such as 1000GP) to increase sample size and variant coverage for imputation of GWAS cohorts. This is already a ubiquitous problem, and there are multiple ways to integrate WGS data sets that require different levels of data sharing and computing power. In this work, we suggest a simple approach that should be feasible for most groups that have sufficient computational resources for GWAS imputation. Our approach is to take two-phased reference panels and reciprocally impute them up to the union set of variants, then use this combined panel for GWAS imputation; we have implemented this functionality in IMPUTE2 (ref. 1) (details are shown in Supplementary Table 1 and Supplementary Note 1).
To evaluate this new functionality, we used a combined 1000GP+UK10K panel to perform imputation with pseudo-GWAS data sets drawn from the UK and Italy (details below). In each of these comparisons, we imputed all available reference variants and stratified them by expected r2, which is a confidence metric produced by IMPUTE2 (also known as ‘info' in the software output). Unlike the true r2 metric, which is usually calculated by masking and imputing ‘truth' genotypes, the expected r2 metric allows direct comparisons of reference panel performance across study populations that have substantially different sets of genotyped truth variants. We have found that predicted r2 values tend to be larger than true r2 values for low-frequency variants (for example, only ~2/3 of variants with expected r2≥0.4 and MAF<5% have true r2≥0.4), so the absolute numbers of high-confidence imputed variants reported in this section should be treated as upper bounds; the emphasis is on qualitative patterns between reference panels and between study populations.
Figure 1b shows how a combined 1000GP+UK10K panel (red) produced by this method performed against each panel separately (1000GP, black; UK10K, blue) when imputing a pseudo-GWAS of UK ancestry. For these evaluations, we used UK10K and 1000GP haplotype panels rephased using SHAPEIT v2, which were previously shown to yield more accurate imputation compared with the corresponding ‘original' haplotypes. The combined and UK10K panels produced very similar numbers of high-confidence (expected r2>0.8) variants at MAFs of 0.5% and higher, implying that the combined panel is neither helpful nor harmful for imputing common and low-frequency variants when a large, population-specific panel is available. On chromosome 20, the combined panel added 2,356 high-confidence rare variants that were not captured by the UK10K panel (MAF<0.5%; 4% increase), which could reflect mutations that have drifted to very low frequencies in the UK but persist on the same haplotype background elsewhere in Europe5,12.
Figure 1c provides the results of a similar evaluation carried out in a population in northern Italy (INCIPE cohort), also based on chromosome 20. The INCIPE cohort was newly genotyped in this study, using Illumina HumanCoreExome-12v1-1 arrays. After stringent QC (Online Supplementary Methods), chromosome 20 genotypes from 6,300 SNPs in 2,145 participants were used to drive imputation with each reference panel. In this data set the UK10K reference panel outperformed the 1000GP panel in all frequency bins, despite the fact that the 1000GP includes a panel (TSI, or ‘Toscani in Italia') that is genetically more similar to the study population. This confirms previous findings13 that reference sample size is often more important than population matching. As before, the combined 1000GP+UK10K panel yielded a larger number of high-confidence imputed variants than the UK10K panel alone—here, the combined panel added 7,466 well-imputed variants with MAF<0.5%, for a 40% increase in rare variants over the UK10K panel (Fig. 1b). These results suggest that it can be especially useful to combine the strengths of multiple panels when a large, population-specific reference set is not available for a particular GWAS population.
Imputation metrics for choosing reference haplotypes
In the course of our analyses, we noticed that some rare variants were imputed well when using the entire UK10K reference panel to drive imputation, yet poorly when using IMPUTE2's khap approximation (all of the results described above are based on using the full reference panel). This approximation reduces the computational cost of imputation by using a region-wide (for example, across a 3MB imputation chunk) Hamming distance metric to reduce the number of reference haplotypes used by a given GWAS haplotype (see also Supplementary Fig. 4). Our investigation of these variants led us to develop a new approximation that uses local (rather than region-wide) haplotype sharing to choose a subset of reference haplotypes (see Supplementary Note 2 for details). This approximation delivers a speed boost similar to that of the existing khap approximation, but it does not sacrifice imputation accuracy at rare and low-frequency variants. For example, Fig. 1d shows the results of imputing the INCIPE pseudo-GWAS data with the UK10K reference panel (see also Supplementary Fig. 5). The full UK10K panel produced the highest accuracy (solid blue line), whereas the khap approximation based on Hamming distance (solid orange line) was less accurate for SNPs with MAF<5%. By contrast, our new approximation based on haplotype tract sharing (dashed orange line) was nearly as accurate as the full reference panel, at ~10% of the computing time (see also Supplementary Fig. 6). All of these strategies for choosing reference haplotypes improved slightly (1–5% increase in mean r2) when the 1000GP haplotypes were added to the UK10K panel, but their relative accuracies remained similar to those shown in Fig. 1d. Further speed improvements are possible for a modest price in accuracy (see Supplementary Note 2).
Discussion
As WGS becomes a standard tool for population and disease genetics, there will be many questions about how to design sequencing studies, how to process the data, how to combine data across studies, and how to limit the computational costs of downstream analysis. With data from one of the most ambitious population sequencing studies to date, we have demonstrated the value of a large, UK-specific reference panel for imputation in British cohorts and in other European populations. Our results show that state-of-the-art phasing methods like SHAPEIT v2 are essential for creating high-quality haplotype panels. Combining WGS data across studies is a desirable goal, and we have implemented an approach in IMPUTE2 that can integrate sets of phased haplotypes to produce a unified reference panel; other strategies for combining WGS data may improve haplotype quality, but our approach has the advantage of being relatively simple and fast. Finally, we have proposed a new approximation that will help reduce the trade-off between imputation speed and accuracy as reference panels continue to grow. The novel strategies we have presented will inform other investigators who wish to use WGS reference panels for imputation, and they will spur additional methods development as population sequencing resources proliferate.
Future efforts to combine multiple large low-coverage sequencing datasets into a substantially larger haplotype resource will likely increase imputation performance, especially at variants with frequencies below 0.1%. We generated a combined reference panel with 42.4 million imputable sites, which is much larger than the 26.6 million imputable sites in the UK10K panel or 32.5 million imputable sites in the 1000GP panel. The UK10K WGS haplotypes for 3,781 samples are available for download from the European Genome-phenome Archive (https://www.ebi.ac.uk/ega/) under managed access conditions (http://www.uk10k.org/data_access). The functionality described in this work is available from the IMPUTE2 website (http://mathgen.stats.ox.ac.uk/impute/impute_v2.html) and the SHAPEIT v2 website (https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html).
Methods
Sample collections
The ALSPAC is a long-term health research project. More than 14,000 mothers enrolled during pregnancy in 1991 and 1992, and the health and development of their children has been followed in great detail ever since8. A random sample of 2,040 study participants was selected for WGS. The ALSPAC Genetics Advisory Committee approved the study and all participants gave signed consent to the study.
The Department of Twin Research and Genetic Epidemiology, is the UK's only twin registry of 11,000 identical and non-identical twins between the ages of 16 and 85 years (ref. 14). The database used to study the genetic and environmental aetiology of age-related complex traits and diseases. The St Thomas's Hospital Ethics Committee approved the study and all participants gave signed consent to the study.
Sequence data production
Low-read depth WGS was performed in the TwinsUK and ALSPAC as part of the UK10K project. Methods for the generation of these data are described in detail as follows6:
Low coverage WGS was performed at both the Wellcome Trust Sanger Institute and the Beijing Genomics Institute (BGI). DNA (1–3μg) was sheared to 100–1,000bp using a Covaris E210 or LE220 (Covaris, Woburn, MA, USA). Sheared DNA was size subjected to Illumina paired-end DNA library preparation. Following size selection (300–500bp insert size), DNA libraries were sequenced using the Illumina HiSeq platform as paired-end 100 base reads according to manufacturer's protocol.
Data generated at the Sanger Institute and BGI were aligned to the human reference separately by the respective centres. The BAM files3 produced from these alignments were submitted to the European Genome-phenome Archive. The Vertebrate Resequencing Group at the Sanger Institute then performed further processing.
Sequencing reads that failed QC were removed using the Illumina GA Pipeline, and the rest were aligned to the GRCh37 human reference, specifically the reference used in Phase 1 of the 1000GP (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz). Reads were aligned using BWA (v0.5.9-r16) (ref. 4). This involved the following steps:
1. Index the reference fasta file:
bwa index -a bwtsw <reference_fasta>
2. For each fastq file:
bwa aln -q 15 -f <sai_file> <reference_fasta> <fastq_file>
3. Create SAM files [sam] using bwa sampe for paired-end reads:
bwa sampe -f <sam_file> <reference_fasta> <sai_files> <fastq_files>
4. Create sorted BAM from SAM. For alignments created at the Sanger this was done using Picard (v1.36; http://picard.sourceforge.net/) SamFormatConverter and samtools (v0.1.11) sort. For alignments created at the BGI, this was done using samtools (v0.1.8) view and samtools sort.
5. PCR duplicates reads in the Sanger alignments were marked as duplicate using the Picard MarkDuplicates, whereas in the BGI alignments they were removed using samtools rmdup.
Further processing to improve SNP and INDEL calling, including realignment around known INDELs, base quality score recalibration, addition of BAQ tags, merging and duplicate marking follows that used for Illumina low coverage data in Phase 1 of the 1000GP5. Software versions used for UK10K for the steps described in that section were GATK version 1.1-5-g6f43284, Picard version 1.64 and samtools version 0.1.16.
SNP and INDEL calls were made using samtools/bcftools (version 0.1.18-r579: https://github.com/samtools/samtools/commit/70c740facc966321754c6bfcc6d61ea056480638)6 by pooling the alignments from 3,910 individual low coverage BAM files. All-samples and all-sites genotype likelihood files (bcf) were created with the samtools mpileup command
samtools mpileup -EDVSp -C50 -m3 -F0.2 -d 8000 -P ILLUMINA -g
with the flags:
C=Coefficient for downgrading mapping quality for reads containing excessive mismatches.
d=At a position, read maximally d reads per input BAM
Variants were then called using the following bcftools command to produce a VCF file7
bcftools view -m 0.9 -vcgN.
For calling on chromosome X and Y, the following settings were applied. The pseudo-autosomal region (PAR) was masked on chromosome Y in the reference fasta file. Male samples were called as diploid in the PAR on chromosome X, and haploid otherwise. No calls were made on chromosome Y for female samples. Diploid/haploid calls were made using the -s option in bcftools view. The PAR regions were: X-PAR1 (60,001-2,699,520); X-PAR2 (154,931,044-155,260,560); Y-PAR1 (10,001-2,649,520); Y-PAR2 (59,034,050-59,363,566). The pipeline (run-mpileup) used to create the calls is available from https://github.com/VertebrateResequencing/vr-codebase/tree/develop.
The observation of spikes in the insertion/deletion ratio in sequencing cycles of a subset of the sequencing runs were linked to the appearance of bubbles in the flow cell during sequencing. To counteract this, the following post-calling filtering was applied. The bamcheck utility from the samtools package was used to create a distribution of INDELs per sequencing cycle. Lanes with INDELs predominantly clustered at certain read cycles were marked as problematic, specifically where the highest peak was 5x bigger than the median of the distribution. The list of problematic lanes included 159 samples. In the next step we checked mapped positions of the affected reads to see if they overlapped with called INDELs, which they did for 1,694,630 called sites. The genotypes and genotype likelihoods of affected samples were then set to the reference genotype unless there was a support for the indel also in a different, unaffected lane from the same sample. In total, 140,163 genotypes were set back to reference and 135,647 sites were excluded by this procedure. Note that this step was carried out on raw, unfiltered calls prior to Variant Quality Score Recalibration (VQSR) filtering.
VQSR8 was used to filter sites. For SNPs, the GATK (version 1.3-21) UnifiedGenotyper was used to recall the sites/alleles discovered by samtools in order to generate annotations to be used for recalibration. Recalibration for the INDELs used annotations derived from the built-in samtools annotations. The GATK VariantRecalibrator was then used to model the variants, followed by GATK ApplyRecalibration, which assigns VQSLOD (variant quality score log odds ratio) values to the variants. For more detailed information on VQSR, see http://www.broadinstitute.org/gsa/wiki/index.php/Variant_quality_score_recalibration. SNPs and INDELs were modeled separately, with parameters given below:
Annotations
SNPs: QD, DP, FS, MQ, HaplotypeScore, MQRankSum, ReadPosRankSum, InbreedingCoeff
INDELSs: MSD, MDV, MSQ, ICF, DP, SB, VDB
Training set
SNPs: HapMap 3.3: hapmap_3.3.b37.sites.vcf, Omni 2.5M chip: 1000G_omni2.5.b37.sites.vcf
INDELs: Mills-Devine, 1000 Genomes Phase I
Truth Set
SNPs: HapMap 3.3: hapmap_3.3.b37.sites.vcf
INDELS: Mills-Devine
Known Set
SNPs: dbSNP build 132: dbsnp_132.b37.vcf
INDELs: Mills-Devine
The truth-set sites are defined as truly showing variation from the reference. VQSLOD scores are calibrated by how many of the truth sites are retained when sites with a VQSLOD score below a given threshold are filtered out. For single-nucleotide variants sites a truth sensitivity of 99.5%, which corresponded to a minimum VQSLOD score of −0.6804 was selected, that is, for this threshold 99.5% of truth sites were retained. For INDEL sites a truth sensitivity of 97%, which corresponded to a minimum VQSLOD score of 0.5939 was chosen. Finally, we also introduced the filter P<10−6 to remove sites that failed the Hardy–Weinberg equilibrium.
The VQSLOD score and other annotations from GATK (BaseQRankSum, Dels, FS, HRun, HaplotypeScore, InbreedingCoeff, MQ0, MQRankSum, QD, ReadPosRankSum, culprit) were copied back to the original samtools calls, excluding annotations which already existed in or did not apply to the samtools VCFs (DP and MQ, AC, AN). Each VCF further contained the filters LowQual (a low-quality variant according to GATK) and MinVQSLOD (Variant's VQSLOD score is less than the cutoff). All sites that did not fail these filters were marked as PASS and brought forward to the genotype refinement stage.
Of the 4,030 samples (1,990 TwinsUK and 2,040 ALSPAC) that were submitted for sequencing, 3,910 samples (1,934 TwinsUK and 1,976 ALSPAC) were sequenced and went through the variant calling procedure. Low-quality samples were identified before the genotype refinement by comparing the samples to their GWAS genotypes using about 20,000 sites on chromosome 20. Comparing the raw genotype calls to existing GWAS data, we removed a total of 112 samples (64 TwinsUK and 48 ALSPAC) because of one or more of the following causes: (i) high overall discordance to SNP array data (>3%; 55 TwinsUK and 36 ALSPAC), (ii) heterozygosity rate>3SD from population mean (1 TwinsUK and 1 ALSPAC), suggesting contamination (iii) no SNP array data available for that sample (7 TwinsUK and 0 ALSPAC) and (iv) sample below 4x mean coverage (1 TwinsUK and 11 ALSPAC). Overall, 3,798 samples (1,870 TwinsUK and 1,928 ALSPAC) were brought forward to the genotype refinement step.
The missing- and low-confidence genotypes in the filtered VCFs were filled out through an imputation procedure with BEAGLE 4 (rev909) (ref. 9).
Additional sample-level QC steps were carried out on refined genotypes, leading to the exclusion of additional 17 samples (16 TwinsUK and 1 ALSPAC) because of one or more of the following causes: (i) non-reference discordance (NRD) with GWAS SNP data>5% (12 TwinsUK and 1 ALSPAC), (ii) contamination identified by multiple relations to other samples (13 TwinsUK and 1 ALSPAC), (iii) failed sex check (3 TwinsUK and 0 ALSPAC). To identify contamination we pruned the WGS data to a set of independent SNPs and calculated genome-wide average identity by state between each pair of samples across the two cohorts. Samples were removed if they had >25 relations with IBS>0.125 (a high number of relationships may indicate contamination). The resulting set of contaminated samples corresponded almost completely to the set of samples with NRD>5%. This left a final set of 3,781 samples (1,854 TwinsUK and 1,927 ALSPAC). These VCF files were submitted to the EGA.
Evaluation of imputation accuracy in the UK10K project
The UK10K final release WGS data of 3,781 samples and 45,492,035 sites was used for creation of haplotype reference WGS data sets. For each chromosome, a summary file was first generated and merged with that of the 1000GP WGS data to identify multi-allelic sites, sites with inconsistent alleles with that of the 1000GP data, and singletons not existing in 1000GP. These sites were excluded to create a new set of VCF files, leaving 26,032,603 sites. The VCF-QUERY tool was used to convert the new VCF files into phased haplotypes and legend files for IMPUTE2. VCF files were converted to binary ped (bed) format and multi-allelic sites excluded, and files were then split into 3MB chunks with ±250kb flanking regions. SHAPEIT v2 was used to re-phrase the haplotypes. Phasing information from the SHAPEIT output was copied back to the original VCF files, with the phase removed for sites missing due to the MAF cutoff. The phased chunks were then recombined with vcf-phased-join from the vcftools package15.
The 1000GP Phase I integrated variant set release (v3) for low-coverage whole-genomes in NCBI build 37 (hg19) coordinates was downloaded from 1000GP FTP site (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20110521/, 23 November 2010 data freezes). This callset includes phased haplotypes for 1,092 individuals and 39,527,072 variants (22 autosome and chromosome X). The haplotypes were inferred from a combination of low-coverage genome sequence data, and they contain SNPs, short INDELs, and large deletions. For each chromosome, a summary file was first generated and merged with that of the UK10K WGS data to identify multi-allelic sites and singletons not polymorphic in UK10K. These sites were excluded to create a new set of VCF files. The final reference panel included all 1,092 samples and 32,449,428 sites. The VCF-QUERY tool was used to convert the new VCF files into phased haplotypes and legend files for IMPUTE2.
A random set of 500 samples passing QC filters was chosen from the TwinsUK (N=1,854) and ALSPAC (N=1,927) WGS data sets. Genotypes for a total of 13,413 sites (corresponding to the content of the Illumina HumanHap610 SNP-array) on chromosome 20 were extracted from the UK10K WGS data in these 1,000 samples.
For the INCIPE study, 6,200 Caucasian participants were randomly chosen from the lists of registered patients of 62 randomly selected general practitioners based in four geographical areas in the Veneto region, North-eastern Italy16. A total of 2,258 samples were genotyped with the HumanCoreExome-12v1-1 platform. A total of 542,585 variants were called using Illumina GENCALL algorithm, 244,594 of which are exonic variants. We conducted further QC evaluation as follows to determine sample and SNP quality. At sample level, we applied the following criteria (i) sample identity was validated through genotyping with an independent typing platform (Sequenom). No samples failed this step. (ii) Twelve pairs of duplicate samples, defined as pairs of individuals with ≥98% concordance genome-wide, were identified. The sample with the lowest call rate of the pair was excluded. (iii) Supplied gender was compared with the genotype-inferred gender (heterozygosity on sample chrX, or is the proportion of chrX SNPs called AB). A Gaussian mixture model was used to find adaptive thresholds Mmax and Fmin (respectively, the maximum male and minimum female heterozygosity on chrX). Overall, 55 samples had chrX heterozygosities that were between Mmax and Fmin, and were excluded from analysis. (iv) Call rate: 90 samples with call rates below 95% were excluded from analysis. (v) 88 samples with autosomal heterozygosity (that is, the proportion of all SNPs with an heterozygous call) score≥3 standard deviations away from the mean were excluded. (vi) Finally, five samples were recommended for exclusion where the normalised magnitude of intensity signal in both channels falls below 0.9. Overall, of the total of 2,258 samples genotyped, 2,145 passed QC filters while 113 samples failed QC filters as indicated above, with some samples failing multiple QC filters. At SNP level, we excluded variants with missingness rate ≥3% or Hardy-Weinberg disequilibrium P<1 × 10−5. We also checked all alleles to confirm that they are on the positive strand of the human genome by comparing alleles against the 1000G and UK10K data. At the end, there were a total of 346,941 polymorphic variants on autosomes, and 8,822 of those on chromosome 20 were retained for analysis.
Pseudo-GWAS panel: for our imputation evaluation, we used 6,300 SNPs on chromosome 20 to mimic a SNP chip in a pseudo-GWAS data set. Before imputation, the two pseudo-GWAS data sets were pre-phased using SHAPEIT v2 (ref. 11) to increase phasing accuracy. The UK10K panel was phased jointly with the entire WGS data set. The INCIPE pseudo-GWAS of 2,145 participants was pre-phased separately.
SHAPEIT v2 was also used for re-phasing the reference haplotypes provided 1000GP and UK10K projects. Per the recommendation of the software, the mean size of the windows in which conditioning haplotypes are defined is set to 0.5MB, instead of 2MB used for pre-phasing GWAS. Owing to the significantly higher number of variants in the WGS data, the re-phasing was conducted by 3MB chunk with 250kb buffering regions, rather than by whole chromosomes as for the pseudo-GWAS. Imputation was carried out on the same chunks with the same flanking regions.
The following three steps were used to merge two WGS reference panels using IMPUTE2 (version 2.3 and later):
Impute the variants that are specific to panel 1 (1000GP) into panel 2 (UK10K).
Impute the variants that are specific to panel 2 (UK10K) into panel 1 (1000GP).
Treat the imputed haplotypes in both panels (with the union of variants from both) as known (that is, take the best-guess haplotypes) and impute the GWAS cohort in the usual way.
The commands for combining haplotypes with the 1000GP are given in Supplementary Note 3.
Imputation of genotypes from the three phased reference panels (UK10K, 1000GP and UK10K+1000GP) into the two test panels was carried out on chromosome 20 split in 3MB chunks with 250kb buffer regions. Imputation was performed using standard parameters with IMPUTE2, for example:
./impute2 \
-m genetic_map_chr20_combined_b37.txt \
-h chr20.uk10k.hap.gz \
-l chr20.uk10k.legend.gz \
-known_haps_g chr20.incipe2gwas.known_haps.gz \
-k_hap 10000 \
-int 3e6 6e6 \
-Ne 20000 \
-buffer 250 \
-use_prephased_g \
-o_gz \
-o chr20.01.incipe2gwas.uk10kRef.impute2
In Fig. 1a,d, the accuracy of imputed variants was calculated as the Pearson correlation coefficient (r2) between imputed genotype dosages in (0–2) and masked sequence genotypes in (0,1,2). The results were stratified into non-overlapping MAF bins for plotting. In Fig. 1b,c, the numbers of variants in different imputation accuracy bins were estimated via the expected r2 (‘info') metric produced by IMPUTE2 (ref. 13). As discussed in the main text, this metric is biased upward relative to the true r2, so the numbers of high-confidence variants in these figures should be interpreted as upper bounds.
Additional information
Accession Codes: UK10K reference haplotypes are available from the European Genome-phenome archive under the accession codes EGAS00001000713 (EGA study) and EGAD00001000776 (EGA dataset) under managed access conditions (http://www.uk10k.org/data_access).
How to cite this article: Huang, J. et al. Improved imputation of low frequency and rare variants using the UK10K haplotype reference panel. Nat. Commun. 6:8111 10.1038/ncomms9111 (2015).
Supplementary Material
Supplementary Figures 1-6, Supplementary Table 1, Supplementary Notes 1-3, and Supplementary References
Acknowledgments
This study makes use of data generated by the UK10K Consortium. The Wellcome Trust provided funding for UK10K (award WT091310). N.S. is supported by the Wellcome Trust (Grant Codes WT098051 and WT091310), the European Commission (EUFP7 EPIGENESYS Grant Code 257082 and BLUEPRINT Grant Code HEALTH-F5-2011-282510) and the NIHR. J.B.R. is funded by the CIHR, CQDM, FRSQ and the Jewish General Hospital. H.-F.Z. is supported by Canadian Institutes of Health Research. J.M. acknowledges support from the ERC (Grant no. 617306). NJT works in a unit supported by the UK Medical Research Council (MRC) (MC_UU_12013/3). For ALSPAC, we are extremely grateful to all the families who took part in this study, the midwives for their help in recruiting them and the whole ALSPAC team, which includes interviewers, computer and laboratory technicians, clerical workers, research scientists, volunteers, managers, receptionists and nurses. For further details concerning the generation of UK10K sequence data, please see extended methods in reference 6. Key contributing studies from the University College London-London School of Hygiene and Tropical Medicine-Edinburgh-Bristol (UCLEB) Consortium appear in the lipids meta-analysis group, but we also acknowledge the support and contribution of the resource more generally. A full description of the collection and its component studies can be found here: DOI: 10.1371/journal.pone.0071345.
Footnotes
Author contributions N.S. and J.M. designed the study. N.S., N.J.T., J.M., J.H. and B.H. wrote the manuscript. B.H. and J.M. developed the software. J.H. and B.H. performed the analyses. Y.M., K.W., J.L.M., P.D., H-F.Z., S.M., J.B.R., R.D. N.J.T., E.T., G.M. and G.G. contributed materials. All authors reviewed and approved the manuscript.
References
- Howie B. N., Donnelly P. & Marchini J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet. 5, e1000529 (2009). [Europe PMC free article] [Abstract] [Google Scholar]
- Frazer K. A. et al.. A second generation human haplotype map of over 3.1 million SNPs. Nature 449, 851–861 (2007). [Europe PMC free article] [Abstract] [Google Scholar]
- Abecasis G. R. et al.. An integrated map of genetic variation from 1,092 human genomes. Nature 491, 56–65 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
- Howie B., Fuchsberger C., Stephens M., Marchini J. & Abecasis G. R. Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. Nat. Genet. 44, 955–959 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
- Howie B., Marchini J. & Stephens M. Genotype imputation with thousands of genomes. G3 (Bethesda) 1, 457–470 (2011). [Europe PMC free article] [Abstract] [Google Scholar]
- The UK10K Consortium. The UK10K project identifies rare variants in health and disease. Nature 10.1038/nature14962 (2015). [Abstract] [Google Scholar]
- Moayyeri A., Hammond C. J., Valdes A. M. & Spector T. D. Cohort profile: twinsUK and healthy ageing twin study. Int. J. Epidemiol. 42, 76–85 (2013). [Europe PMC free article] [Abstract] [Google Scholar]
- Golding J., Pembrey M. & Jones R. ALSPAC–the avon longitudinal study of parents and children. I. Study methodology. Paediatr. Perinat. Epidemiol. 15, 74–87 (2001). [Abstract] [Google Scholar]
- Li H. et al.. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009). [Europe PMC free article] [Abstract] [Google Scholar]
- Browning B. L. & Browning S. R. A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. Am. J. Hum. Genet. 84, 210–223 (2009). [Europe PMC free article] [Abstract] [Google Scholar]
- Delaneau O., Zagury J. F. & Marchini J. Improved whole-chromosome phasing for disease and population genetic studies. Nat. Methods 10, 5–6 (2013). [Abstract] [Google Scholar]
- Jewett E. M., Zawistowski M., Rosenberg N. A. & Zollner S. A coalescent model for genotype imputation. Genetics 191, 1239–1255 (2012). [Europe PMC free article] [Abstract] [Google Scholar]
- Marchini J. & Howie B. Genotype imputation for genome-wide association studies. Nat. Rev. Genet. 11, 499–511 (2010). [Abstract] [Google Scholar]
- Moayyeri A., Hammond C. J., Hart D. J. & Spector T. D. The UK adult twin registry (TwinsUK Resource). Twin Res. Hum. Genet. 16, 144–149 (2013). [Europe PMC free article] [Abstract] [Google Scholar]
- Danecek P. et al.. The variant call format and VCFtools. Bioinformatics. 27, 2156–2158 (2011). [Europe PMC free article] [Abstract] [Google Scholar]
- Gambaro G. et al.. Prevalence of CKD in northeastern Italy: results of the INCIPE study and comparison with NHANES. Clin. J. Am. Soc. Nephrol. 5, 1946–1953 (2010). [Europe PMC free article] [Abstract] [Google Scholar]
Full text links
Read article at publisher's site: https://doi.org/10.1038/ncomms9111
Read article for free, from open access legal sources, via Unpaywall:
https://www.nature.com/articles/ncomms9111.pdf
Citations & impact
Impact metrics
Citations of article over time
Alternative metrics
Article citations
Genetic basis of right and left ventricular heart shape.
Nat Commun, 15(1):9437, 14 Nov 2024
Cited by: 0 articles | PMID: 39543113 | PMCID: PMC11564811
Multiomic integration analysis identifies atherogenic metabolites mediating between novel immune genes and cardiovascular risk.
Genome Med, 16(1):122, 24 Oct 2024
Cited by: 0 articles | PMID: 39449064 | PMCID: PMC11515386
Population-specific putative causal variants shape quantitative traits.
Nat Genet, 56(10):2027-2035, 03 Oct 2024
Cited by: 1 article | PMID: 39363016 | PMCID: PMC11525193
Mapping biological influences on the human plasma proteome beyond the genome.
Nat Metab, 6(10):2010-2023, 26 Sep 2024
Cited by: 1 article | PMID: 39327534 | PMCID: PMC11496106
Integration of estimated regional gene expression with neuroimaging and clinical phenotypes at biobank scale.
PLoS Biol, 22(9):e3002782, 13 Sep 2024
Cited by: 0 articles | PMID: 39269986 | PMCID: PMC11424006
Go to all (196) article citations
Data
Data behind the article
This data has been text mined from the article, or deposited into data resources.
BioStudies: supplemental material and supporting data
European Genome-Phenome Archive (2)
- (1 citation) European Genome-Phenome Archive - EGAD00001000776
- (1 citation) European Genome-Phenome Archive - EGAS00001000713
Similar Articles
To arrive at the top five similar articles we use a word-weighted algorithm to compare words from the Title and Abstract of each citation.
Rare variant genotype imputation with thousands of study-specific whole-genome sequences: implications for cost-effective study designs.
Eur J Hum Genet, 23(7):975-983, 08 Oct 2014
Cited by: 69 articles | PMID: 25293720 | PMCID: PMC4463504
Genotype imputation performance of three reference panels using African ancestry individuals.
Hum Genet, 137(4):281-292, 10 Apr 2018
Cited by: 26 articles | PMID: 29637265 | PMCID: PMC6209094
Improved imputation quality of low-frequency and rare variants in European samples using the 'Genome of The Netherlands'.
Eur J Hum Genet, 22(11):1321-1326, 04 Jun 2014
Cited by: 68 articles | PMID: 24896149 | PMCID: PMC4200431
Genotype Imputation in Genome-Wide Association Studies.
Curr Protoc Hum Genet, 102(1):e84, 01 Jun 2019
Cited by: 17 articles | PMID: 31216114
Review
Funding
Funders who supported this work.
British Heart Foundation (3)
Tbx1 and cardiovascular morphogenesis: genetic networks and tissue interactions (renewal: years 11-15)
Peter Scambler, University College London
Grant ID: RG/10/13/28570
High resolution analysis of biological and environmental determinants of cardiovascular risk factors and cardiovascular events in older women
Professor Juan Casas, London School of Hygiene and Tropical Medicine
Grant ID: PG/13/66/30442
Signalling pathways in cardiogenesis
Shoumo Bhattacharya, University of Oxford
Grant ID: RG/10/17/28553
CIHR
Economic and Social Research Council (1)
Grant ID: ES/M001660/1
European Research Council (1)
ADVANCED STATISTICAL METHODS FOR HIGH-DIMENSIONAL GENETIC STUDIES (HIGEN)
Prof Jonathan Marchini, University of Oxford
Grant ID: 617306
Medical Research Council (9)
The Avon Longitudinal Study of Parents and Children: An international resource for population genomics and lifecourse epidemiology. Core Programme Support 2011-2015 and Core programme support 2014-2019
Professor Nicholas Timpson, The Wellcome Trust Ltd
Grant ID: MC_PC_15018
Causal analyses, statistical efficiency and phenotypic precision through study design: Genotype based sampling
Professor Nicholas Timpson, University of Bristol
Grant ID: MC_UU_12013/3
MRC Centre for Neuropsychiatric Genetics and Genomics
Professor Sir Michael Owen, Cardiff University
Grant ID: MR/L010305/1
Enabling technologies
Professor Giles Yeo, University of Cambridge
Grant ID: MC_UU_12012/5
Molecular Genetics of Schizophrenia
Professor Sir Michael Owen, Cardiff University
Grant ID: G0800509
Aetiology of type 2 diabetes and related metabolic disorders
Professor Nicholas Wareham, University of Cambridge
Grant ID: MC_UU_12015/1
From Mendelian Randomization to Hypothesis Free Causal Inference
Professor George Davey Smith, University of Bristol
Grant ID: MC_UU_12013/1
Using Genetics to Identify Causal Pathways that Influence Bone Related Phenotypes in Children and Young Adults.
Professor David Evans, University of Bristol
Grant ID: MC_UU_12013/4
Grant ID: MC_UU_12012/5/B
National Institute for Health Research (NIHR) (8)
Grant ID: NF-SI-0510-10268
Grant ID: NF-SI-0508-10198
Grant ID: NF-SI-0513-10109
Grant ID: NF-SI-0513-10008
Head injury: monitoring and optimising cerebral metabolism to improve outcome
Professor Peter Hutchinson, University of Cambridge
Grant ID: NIHR-RP-R3-12-013
Grant ID: NF-SI-0507-10380
Grant ID: NF-SI-0514-10027
Grant ID: NF-SI-0514-10176
Wellcome Trust (12)
Disorders of nuclear hormone synthesis & action: genetics and pathophysiology.
Prof Krishna Chatterjee, University of Cambridge
Grant ID: 095564
Molecular and neural basis of obesity
Professor Ismaa Sadaf Farooqi, University of Cambridge
Grant ID: 098497
Grant ID: WT091310
Grant ID: WT098051
Genetic dissection of mechanisms linking cell dysfunction, insulin resistance & major human disease
Prof Robert Semple, University of Cambridge
Grant ID: 098498
LIPODYSTROPHY-A PARADIGM FOR ELUCIDATING PATHOGENIC MECHANISMS IN THE METABOLIC SYNDROME.
Prof David Savage, University of Cambridge
Grant ID: 091551
Cellular Medicine at the Cambridge Institute for Medical Research.
Prof Gillian Griffiths, University of Cambridge
Grant ID: 100140
Insulin Resistance: Lessons from Extreme Phenotypes.
Prof Sir Stephen O'Rahilly, University of Cambridge
Grant ID: 095515
Embedding Biobanks as Tools for Genomic Medicine.
Dr Jane Kaye, University of Oxford
Grant ID: 096599
The Avon Longitudinal Study for Parents and Children: An international resource for population genomic and lifecourse epidemiology. Core programme support 2014-2019.
Professor George Davey Smith, University of Bristol
Grant ID: 102215
Institute of Metabolic Science.
Prof Sir Stephen O'Rahilly, University of Cambridge
Grant ID: 100574
Grant ID: 104036/Z/14/Z
