Evaluating the quality of the 1000 genomes project data

Background

Data from the 1000 Genomes project is quite often used as a reference for human genomic analysis. However, its accuracy needs to be assessed to understand the quality of predictions made using this reference. We present here an assessment of the genotyping, phasing, and imputation accuracy data in the 1000 Genomes project. We compare the phased haplotype calls from the 1000 Genomes project to experimentally phased haplotypes for 28 of the same individuals sequenced using the 10X Genomics platform.

Results

We observe that phasing and imputation for rare variants are unreliable, which likely reflects the limited sample size of the 1000 Genomes project data. Further, it appears that using a population specific reference panel does not improve the accuracy of imputation over using the entire 1000 Genomes data set as a reference panel. We also note that the error rates and trends depend on the choice of definition of error, and hence any error reporting needs to take these definitions into account.

Conclusions

The quality of the 1000 Genomes data needs to be considered while using this database for further studies. This work presents an analysis that can be used for these assessments.

The 1000 Genomes Project (1000GP) was designed to provide a comprehensive description of human genetic variation through sequencing multiple individuals [1,2,3]. Specifically, the 1000GP provides a list of variants and haplotypes that can be used for evolutionary, functional and biomedical studies of human genetics. Over the three phases of the 1000GP, a total of 2504 individuals across 26 populations were sequenced. These populations were classified into 5 major continental groups: Africa (AFR), America (AMR), Europe (EUR), East Asia (EAS), and South Asia (SAS). The 1000GP data was generated using a combination of multiple sequencing approaches, including low coverage whole genome sequencing with mean depth of 7.4X, deep exome sequencing with a mean depth of 65.7X, and dense microarray genotyping. These sequences were used for calling genotypes and generating the variant calls. In addition, a subset of individuals (427) including mother-father-child trios and parent-child duos were deep sequenced using the Complete Genomics platform at a high coverage mean depth of 47X. The project involved characterization of biallelic and multiallelic SNPs, indels, and structural variants.

Given the low depth of (sequencing) coverage for most 1000GP samples, it is unclear how accurate the imputed haplotypes are, especially for rare variants. We quantify this accuracy directly by comparing imputed genotypes and haplotypes based on low-coverage whole-genome sequence data from the 1000GP with highly accurate, experimentally determined haplotypes from 28 of the same samples. Additional motivation for our study is given below.

Phasing

It is important to understand phase information in analyzing human genomic data. Phasing involves resolving haplotypes for sites across individual whole genome sequences. The term ‘diplomics’ [4] has been coined to describe “scientific investigations that leverage phase information in order to understand how molecular and clinical phenotypes are influenced by unique diplotypes”. The diplotype shows effects in function and disease related phenotypes. Multiple phenomena like allele-specific expression, compound heterozygosity, inferring human demographic history, and resolving structural variants requires an understanding of the phase of available genomic data. Phased haplotypes are also required as an intermediate step for genotype imputation.

Phasing methods can be categorized into methods which use information from multiple individuals and those which rely on information from a single individual [5]. The former are primarily computational methods, while the latter are mostly experimental approaches. Some computational approaches use information from existing population genomic databases and can be used for phasing multiple individuals. These, however, may be unable to correctly phase rare and private variants, which are not represented in the reference database used. On the other hand, some methods use information from parents or closely related individuals. These have the advantage of being able to use Identical-By-Descent (IBD) information, and allow long range phasing, but require sequencing of more individuals, which adds to the cost. Commonly used computational phasing methods are: BEAGLE [6], SHAPEIT [7, 8], EAGLE [9, 10] and IMPUTE v2 [11].

Experimental phasing methods, on the other hand, often involve separation of entire chromosomes followed by sequencing of short segments, which can then be computationally reconstructed to generate entire haplotypes. These methods do not need information from individuals other than the one being sequenced. They involve genotyping being performed separately from phasing. These methods fall into two broad categories, namely dense and sparse methods [12]. Dense methods resolve haplotypes in small blocks in great detail, where all variants in a specific region are phased. However, they do not inform the phase relationship between the haplotype blocks. These involve diluting high molecular weight DNA fragments such that fragments from at most one haplotype are present in each unit. Sparse methods can resolve phase relationships across large distances, but may not inform on the phase of each variant in a chromosome. In these methods, a low number of whole chromosomes is compartmentalized such that only one of each pair of haplotypes is present in each compartment. These compartmentalizations are followed by sequencing to generate the haplotypes.

In this work, we use phased haplotypes generated using the 10X Genomics method which uses linked-read sequencing [13]. This method can be best classified as a dense phasing method. Most of the SNPs (~ 99%) are phased. One nanogram of high molecular weight genomic DNA is distributed across 100,000 droplets. This DNA is barcoded and amplified using polymerase. This tagged DNA is released from the droplets and undergoes library preparation. These libraries are processed via Illumina short-read sequencing. A computational algorithm is then used to construct phased haplotypes based on the barcodes. This method has been shown to have the lowest error rate (0.064%) [14]. This error rate is considerably lower than the error rate we observe for the 1000 Genomes phasing (as reported in our Results).

Imputation

Imputation involves the prediction of genotypes not directly assayed in a sample of individuals. Experimentally sequencing genomes to a high coverage is an expensive process. Low coverage sequencing or arrays can be used as low-cost methods for sequencing. However, these methods may lead to uncertainty in estimated genotypes (low coverage sequencing) or missing genotype values for untyped sites (arrays). Imputation can be used to obtain genotype data for missing positions using reference data and known data at a subset of positions in individuals which need to be imputed. Imputation is used to boost the power of GWAS studies [15], fine mapping a particular region of a chromosome [16], or performing meta-analysis [17], which involves combining reference data from multiple reference panels.

Imputation uses a reference panel of known haplotypes with alleles known at a high density of haplotyped positions. A study/inference panel genotyped at a sparse set of positions is used for sequences which need to be imputed. A basic conceptual description of imputation involves phasing genotypes at genotyped positions in the study/inference panel, followed by matching haplotypes which match in the genotyped positions [11]. Various imputation algorithms perform these steps sequentially and iteratively or simultaneously, while others further improve on this basic approach by including probabilistic modeling.

Factors affecting the quality of the phasing and imputation are (1) size of reference panel (2) density of SNPs in reference panel (3) accuracy of called genotypes in the reference panel (4) degree of relatedness between sequences in reference panel and study sequences (5) ethnicity of the study individuals in comparison with the available reference data and (6) allele frequency of the site being phased or imputed [5].

Multiple methods have been developed for genotype imputation [18]. MACH [19, 20], minimac [21], BEAGLE [6], and IMPUTE v2 [11] are some widely used methods for imputation.

An analysis of the imputation accuracy for the HapMap project has been performed about a decade ago [22]. The 1000 Genomes project has performed a similar analysis with the WGS data sequenced with Complete Genomics [3]. We present here a detailed alternative assessment of the quality of phasing and imputation for the 1000 Genomes database comparing with high coverage experimentally phased sequences sequenced using a new method for experimentally resolving haplotypes, particularly as a function of minor allele frequency and inter-SNP distances for biallelic SNPs.

The 1000 Genomes project chromosome-specific VCFs for the GRCh38 assembly contain between 7.07 M (chr2) to 1.1 M (chr22) variants over all the 2504 individuals. After filtering for biallelic SNPs, phased, filtered for PASS, removing indels, we are left with 6.78 M (chr2) to 1.05 M (chr22) variants. The experimentally phased data from the 10X Genomics platform has different numbers of called variants for each sequenced individual. For chromosome 1, the number of called variants varies from 414 K to 494 K across the 28 individuals, while, for chromosome 22, the number of called SNPs varies from 104 K to 120 K. After performing a similar filtering for the experimental data, the number of biallelic PASS phased SNPs ranges between 298 K and 357 K for chromosome 1 and 64 K and 75 K for chromosome 22.

The SNPs from the experimentally phased VCFs (Fig. 1a), averaged over continent groups show that the vast majority of SNPs in this selection have high continent-specific MAF values (> 5%). However, if we look at all the SNPs in the 1000 Genomes Data (filtered for biallelic PASS phased SNPs) as a function of continent-specific MAF, the distribution we observe has a very different trend. There is a significant over-representation of the very low continent-specific MAF SNPs (< 0.1%), ∼ 5 ∗ 10⁷, as compared to all the subsequent higher MAF SNPs, which all range < 1 ∗ 10⁷.

Distribution of SNPs as a function of continent-specific minor allele frequencies a only experimental SNPs b all 1000 Genomes SNPs

Full size image

These discrepancies between the numbers in the 1000 Genomes data and in the experimentally phased data, as well as the differing trends as a function of MAF occur because the 1000 Genomes data includes a SNP if even one individual in the 2504 individuals has a variant (heterozygous or homozygous-alternate) at that position while the experimental data includes a SNP only if that particular individual has a variant (heterozygous or homozygous-alternate) at that position. This results in a much larger number of overall SNPs being present in the 1000 Genomes data as compared to the experimental and also the majority of the 1000 Genomes SNPs having extremely low MAF, as those would occur only in one or a few individuals.

Genotyping error

Genotyping error is computed comparing the 1000 Genomes genotypes with the experimental genotypes. The experimental genotypes for all SNPs not present in the experimental VCF for each individual are assumed to be homozygous reference. Mismatched genotypes are counted as errors. Figure 2a looks at the errors (fraction of genotypes which are incorrect) for the experimental VCF positions as a function of the continent-specific minor allele frequencies. There is higher error at the population invariant sites (MAF = 0.0%) in the African and American populations than the European, East Asian and South Asian populations. This correlates with a lower total number of population invariant SNPs in those continents (Fig. 1a). For non-invariant SNPs, we observe, as expected, a decreasing error rate with increasing minor allele frequency, to a < 2% error genotyping error rate for the SNPs with minor allele frequencies > 1%.

Genotyping error a in the experimental VCF positions (non-hom ref. SNPs) as a function of continent-specific minor allele frequency averaged over all chromosomes over all individuals in each continent b in experimental VCF positions comparing SNPs with homozygous alternate vs heterozygous calls in the experimental data c false positive vs false negative rates (defined in text) for all 1000 Genomes SNPs

Full size image

Within these errors in the experimental SNPs, we observe significantly different rates for SNPs which are heterozygous vs homozygous reference in the experimental data (Fig. 2b). The error rate for SNPs which are homozygous alternate in the experimental data is 1.5 times the error rate for the SNPs which are heterozygous in the experimental data.

Comparing false positive (sites non-homozygous reference in 1000 Genomes data and homozygous reference in the experimental data) vs false negative (sites homozygous reference in 1000 Genomes data and non-homozygous reference in the experimental data) error rates for all 1000 Genomes sites (Fig. 2c), we see that the East Asian and South Asian populations both have mostly low false positive rates, but show a wide range (factor of 2) of false negative rates, while showing only a ~ 15% variation in the false positive rates for most individuals. In contrast, the African individuals mostly have relatively low false negative rates, but have among the highest false positive rates. This indicates that the sequencing in the 1000 Genomes project has over called non-homozygous reference variants in African individuals compared to the rest, and over called SNPs as homozygous reference in some of the East and South Asian individuals.

Phasing

Phasing errors are all analyzed for overall 1000 Genomes minor allele frequencies, not continent specific MAFs. Comparing the switch error across individual chromosomes (Fig. 3), we observe that the switch error ranges between 20 and 30% for the rare MAF (< 0.1%) SNPs, falling to < 5% for SNPs with MAFs 1–5%. The majority of SNPs, which fall in the MAF > 5% category, have an error < 2.5%. However, a comparatively higher switch error at larger MAF values (> 5%) is observed for chromosome 21. This plot (Fig. 3) shows only a subset of chromosomes for a single individual (GM18552), but this trend is observed for all other chromosomes and individuals studied.

Switch error as a function of Minor Allele Frequencies for different individual chromosomes. Chromosome 21 shows higher switch error for large MAF values

Full size image

Figure 4a shows the total switch error for each of the individuals. The total switch errors for all the individuals studied go up to ∼ 2.5%. The switch errors for the East Asian individuals are grouped together, while those for the South Asian individuals show greater variability. This is in line with the general observation that South Asian populations have an overall greater heterogeneity than do East Asian populations, which some of the authors have observed in ongoing studies with hundreds of individuals [J. Wall, Unpublished data].

Switch error a Total switch error (number of switches in experimental SNPs/total number of experimental SNPs) for each individual b Switch error as a function of Minor Allele Frequencies averaged over all individuals in each continent. c Switch error as a function of Minor Allele Frequencies for all individuals colored by continent

Full size image

Analyzing the switch error as a function of minor allele frequency averaged over all chromosomes of all individuals of a population (Fig. 4b), we observe low switch error, < 5%, for low minor allele frequencies (MAF) (1–5%). For rare SNPs with MAF (0.2–1%), the switch error is ∼ 5–10%. For extremely rare minor allele SNPs, i.e. MAF < 0.2%, the error is much higher, i.e. 15–35%. For all higher MAF values (> 5%), the error is < 2.5%. The average error rate for the individuals from the African populations is almost the same over the range of MAF values > 0.1%.

As observed in Fig. 4c, the differences in the error rates between individuals decrease with increasing minor allele frequency. Individuals from South Asia show a larger variation in error as a function of MAF as compared to individuals from East Asia. The individuals from the African populations have the lowest switch error over the range of MAF values. Individual NA20900, an individual from the Gujarati Indians in Houston (GIH) population has the lowest switch error as a function of minor allele frequency for the low MAF SNPs. This individual is not part of a trio in the 1000GP data, and further analysis is required to ascertain why it shows much lower switch error as compared to the other individuals studied. One possible explanation is that the current limited sampling of only 11 individuals from the South Asian population is not capturing the full spread of error rate variation, and including more individuals might show more individuals with comparable low error rates.

We also analyzed phasing error as a function of the distances between SNPs (Fig. 5). The phasing error increases as a function of the inter-SNP distance, i.e. SNPs which are further apart are more likely to be out of phase with each other. The within population trends are the same as for switch error vs MAF, where the individuals from South Asia show a larger spread as compared to the individuals from East Asia. Individual NA20900 shows the lowest error rate, same as for the comparison of error vs MAF (Fig. 4c).

Switch error as a function of inter-SNP distance a Switch error as a function of inter-SNP distances averaged over individuals in each continent. b Switch error as a function of inter-SNP distances for all individuals colored by continent

Full size image

Comparing the switch error as a function of MAF vs. the switch error as a function of inter-SNP distance, we see that the individuals from the African populations show distinctly opposite trends. For low MAF SNPs, the error is the lowest averaging over the African individuals, while across the range of inter-SNP distances, the average over the African individuals was the highest error. The reason this occurs can be understood from the fact that there are a higher number of low MAF SNPs in the African individuals in the experimental data (Fig. 1a), as well as an overall higher number of SNPs in those individuals, leading to a higher SNP density for these individuals. In addition, there is less linkage disequilibrium (LD) in the individuals from the African populations, which would make it harder to phase them accurately [23, 24]. Hence, pairs of SNPs are more likely to be out of phase with each other, leading to higher switch error as a function of inter-SNP distance.

Imputation

Imputation error is computed as the fraction of SNPs with incorrectly imputed genotypes (genotype discordance). However, depending on the subset of SNPs under consideration, the error can be computed in two different ways, (1) fraction of experimental SNPs incorrectly imputed and (2) fraction of all 1000GP SNPs incorrectly imputed. In the case of the second definition of error, the experimental calls for all the positions not in the experimental VCFs are assumed to be homozygous-reference.

Figure 6a shows the total imputation error in the experimental SNPs while Fig. 6b shows the total imputation error in the 1000GP SNPs for each of the individuals. The total imputation errors in the experimental SNPs for all the individuals studied go up to ∼ 4%. For this subset of SNPs, the two American individuals have the among the highest imputation errors. The imputation errors for the East Asian individuals are grouped together, while those for the South Asian individuals show greater variability. This agrees with our observations for the switch error (Fig. 4a). In the 1000GP SNPs, on the other hand, since we are looking at a much larger set of SNPs, most of which are homozygous-reference in any given individual, we see a much smaller error < ∼ 1%.

Total imputation error a Total imputation error in experimental SNPs (number of incorrect genotypes in all experimental SNPs/total number of experimental SNPs) for each individual b Total imputation error in all 1000GP SNPs (number of incorrect genotypes in all 1000GP SNPs/total number of 1000GP SNPs) for each individual

Full size image

Imputation error in experimental SNPs

Figure 7a shows the imputation error rates as function of the continent-specific minor allele frequency. The continent invariant positions (MAF = 0.0%) are imputed almost as accurately as the high MAF (> 5% in 3 populations, and > 1% in two populations) SNPs. In these positions, we make the same observation as we did for the original genotyping in the 1000 genomes reference data (Fig. 2a), i.e. the errors in the European, East Asian and South Asian individuals for these continent invariant positions are lower than those for the American and African individuals. For the very rare SNPs, i.e. MAF < 0.2%, the error is as high as ∼ 60%. These extremely high error rates are only observed in the American individuals and a few of the South Asian individuals. While this error rate seems high, a likely explanation for that is that the imputation method infers each allele by finding the most likely haplotype match from the reference database for the individual being imputed [11]. In the case of a SNP with a rare variant, the best matching haplotypes are likely to contain the reference allele, leading to a prediction of homozygous reference genotype at that position. However, the SNPs in the experimental VCFs only include positions for which there is a non-homozygous reference genotype for that particular individual. As a result, any prediction of homozygous reference genotype is going to be counted as an error, leading to comparatively high error rates at these very low MAF values. For the rest of the individuals, the error rates are < 50%. In the mid-range of MAF values, i.e. 0.2 to 1%, the errors range between 10 and 20%. The SNPs with higher MAF values are fairly accurate, with errors < 2% for common SNPs (MAF > 5%). This can also be seen looking at all the individuals separately (Fig. 7b). The South Asian (Gujarati in Houston, Texas) individual NA20900 still shows the lowest error rate as a function of MAF for imputation, just as it does for the switch error (Fig. 4c). Out of the imputed experimental SNPs, a very small fraction have low imputation INFO scores (Additional file 1: Figure S1a). However, most of those are SNPs which are imputed incorrectly, hence filtering out low INFO score SNPs gives much smaller error rates throughout the range of MAF values (Additional file 1: Figure S2b).

Imputation accuracy experimental VCF positions a Imputation error in the experimental SNPs as a function of Minor Allele Frequencies averaged over individuals in each continent. b Imputation error in the experimental SNPs as a function of Minor Allele Frequencies for all individuals colored by continent

Full size image

Imputation error in all 1000GP SNPs

Computing the error using all the 1000GP SNPs, we see a different trend for the errors as a function of minor allele frequency (Fig. 8a, b). The invariant sites have very low errors ~ 10^− 4. For the variant sites, the errors increase as a function of minor allele frequency, as opposed to decreasing as they do in the experimental only SNPs. The reason this happens is that contrasting the number of experimental SNPs (Fig. 1a) with the numbers of all 1000GP SNPs (Fig. 1b), while the number of low MAF SNPs is 1–2 orders of magnitude less than the number of SNPs with MAF > 5% in the experimental data, the number of very low MAF SNPs is 2–10 times greater than the number of SNPs with MAF > 5% in the whole 1000 Genomes data. The vast majority of the very low MAF SNPs in the whole 1000 Genomes data are homozygous-reference, since those SNPs show variation in only one or very few 1000 Genomes individuals. Hence, imputation predictions get most of those positions correct in most of the individuals. As a result, the fraction of those very rare SNPs which are predicted incorrectly is much lower when considering all the 1000 Genomes SNPs as compared to only considering the experimental SNPs, where most of the SNPs are high MAF SNPs. However, it is important to note that a lot of the low MAF SNPs have low INFO scores for imputation (Additional file 1: Figure S1b). Hence filtering out SNPs with low INFO scores shows a decreasing error rate with increasing MAF, as is expected (Additional file 1: Figure S3b).

Imputation accuracy all 1000GP SNPs a Imputation error in all the 1000 Genomes positions as a function of Minor Allele Frequencies averaged over individuals in each continent. b Imputation error in all the 1000 Genomes positions as a function of Minor Allele Frequencies for all individuals colored by continent

Full size image

Consistent with the observations for the experimental only SNPs, at very rare SNPs (MAF < 0.2%), the American individuals still have the highest error rate. The individuals from the South Asian populations still show a greater spread than those from the East Asian populations. Individual NA20900 still shows the lowest error rate as with previous observations.

An alternative measure of imputation accuracy is genotype r². Figure 9 shows the r² as function of the alternate allele frequency (AAF) (as opposed to minor allele frequencies). This enables comparison to the imputation accuracies reported in the 1000GP phase 3 paper [3], and we see higher accuracies for EAS individuals and lower accuracies for AMR individuals at very low alternate allele frequencies compared to those previously reported values. The accuracies reported for SNPs with AAF > 1% are consistent with the previously reported values in the 1000GP phase 3 paper. Consistent with the observations in genotype discordance, the r² values show the least accuracy for the American individuals at low alternate allele frequencies.

Imputation accuracy all 1000GP SNPs r² for allele frequency bins

Full size image

Comparison of reference panels

Here, we compare the imputation errors resulting from using different reference panels for imputation. A continent-specific reference panel for the individual of interest, a reference panel which includes all of the 1000 Genomes individuals, and a continent-specific reference panel for a different continent from the one from which the individuals are, are chosen. The minor allele frequencies used here are for all the overall 1000 Genomes minor allele frequencies, instead of a continent-specific minor allele frequency, since we want to understand the impact of the choice of reference panel, and continent-specific MAFs would not align with the whole reference or the reference from another continent. In this case, we used the South Asian reference panel as the different continent panel and estimated imputation accuracies for all the other individuals, using a reference panel corresponding to that individual’s continent group, the South Asian reference panel, and the whole 1000G reference panel.

The observed result for experimental only SNPs (Fig. 10a) when comparing reference panels for the AFR, AMR, EUR, and EAS individuals is very similar when looking at all 1000 Genomes SNPs (Fig. 10b). The imputation accuracy when using the entire 1000 Genomes data as a reference panel gives almost identical accuracy as using a continent specific reference panel corresponding to the individuals in 3 of the 4 continent groups. For the AMR individuals, however, there is a marked improvement in using the full 1000G reference panel than the AMR specific reference panel. The error while using an incorrect reference panel, in this case the SAS panel, however, is a factor of 2 or more greater than the error when using the appropriate reference, or when using the whole 1000 Genomes reference panel. In particular, the choice of the SAS panel gives significantly the highest error rate for the AFR individuals. The trend of error as a function of MAF for all 1000G SNPs is, again, the opposite of what was observed when looking at only the experimental SNPs, as discussed previously.

Imputation error as a function of Minor Allele Frequencies for AFR (red), AMR (blue), EUR (black), and EAS (green) individuals comparing the continent specific reference panel (solid lines + circles), a different continent specific panel (SAS, dotted lines + squares), and the entire 1000G reference panel (dashed lines + triangles) a experimental SNPs b All 1000 Genomes SNPs

Full size image

The 1000 Genomes Project data have been widely used as a reference for estimating continent-specific allele frequencies, and as a reference panel for phasing and imputation studies. Since the project’s design involved low-coverage (~7X) sequencing for most of the samples, it was unknown a priori how accurate the 1000GP’s genotype and haplotype calls were, especially for rare variants. This accuracy obviously directly impacts the usefulness of the 1000GP data. While some quantification of imputation accuracy in the 1000GP has been performed before [3], with the advent of inexpensive, commercial platforms for experimentally phasing whole genomes, it is possible to directly quantify the genotype and haplotype error rates of the 1000GP data.

Our comparison of 28 experimentally phased genomes with the 1000GP data found that the latter is highly accurate for common and low-frequency variants (i.e., MAF ≥ 0.01). As expected, accuracy declined with decreasing MAF, with rare variants (MAF < 0.01) not reliably imputed onto haplotypes. Surprisingly though, the genotype calls were reasonably accurate even for rare variants. This observation may not generalize to other low-coverage sequencing studies due to the complicated and labor-intensive protocol used for variant calling in the 1000GP. We conclude that the 1000GP data is best used as a reference panel for imputing variants with MAF ≥ 0.01 into populations closely related to the 1000GP groups, and is probably of limited utility for imputation in rare variant association studies. Larger subsequent imputation panels, such as the one generated by the Haplotype Reference Consortium (HRC) [25], are likely much more useful for imputing rare variants, at least in well-studied European populations. However, even this large reference panel may be of limited usefulness for imputation into other human groups. While our results suggest that using a region-specific reference panel (for the correct region) for imputation is only slightly worse than using a worldwide panel, the choice of an incorrect regional panel makes the imputation considerably worse. So, large European-based haplotype reference panels will be of limited utility for imputing variants into East Asian, South Asian, or African-American genomes, while imputation studies involving understudied groups such as Middle Easterners, Melanesians or Khoisan are likely to have error rates substantially higher than what was observed in our study. This is a consequence of the fact that most rare variants are region-specific; imputation only works when the variant being imputed shows up often enough in the reference panel. In summary, while the 1000GP and HRC provide valuable genomic resources that can augment the power of GWAS in groups with European ancestry, additional large-scale genome sequencing of diverse human populations will be necessary to obtain comparable benefits of imputation in genetic association studies of non-European groups.

Finally, we note that the absolute error rate varied by an order of magnitude, depending on the specific definitions of error that were used. This highlights the importance of definitional clarity in studies that evaluate the accuracy of genomic resources.

Input data

Processed VCFs were downloaded from the 1000 Genomes website. This data is available for each chromosome separately. To obtain agreement with the experimental data, 1000 Genomes VCFs corresponding to the GRCh38 assembly were downloaded. Experimental data was sequenced using the 10X Genomics platform for 28 individuals from the 1000 Genomes project. Thirteen of these individuals were processed at UCSF [26] and sequenced at Novogene, while the remaining individuals were processed and sequenced at Genentech. The populations from which each of the individuals come (as listed in the Coriell Catalog) are:

• South Asia (SAS):

  • ○ Gujarati Indians in Houston, Texas, USA (HapMap) [GIH] - GM21125*, NA20900, NA20902

  • ○ Punjabi in Lahore, Pakistan [PJL] - HG03491, HG03619

  • ○ Sri Lankan Tamil in the UK [STU] - HG03679, HG03752, HG03838*

  • ○ Indian Telugu in the UK [ITU] - HG03968

  • ○ Bengali in Bangladesh [BEB] - HG04153, HG04155

• East Asia (EAS):

  • ○ Han Chinese in Beijing, China (HapMap) [CHB] - GM18552*, NA18570, NA18571

  • ○ Chinese Dai in Xishuangbanna, China [CDX] - HG00851*, HG01802, HG01804

  • ○ Kinh in Ho Chi Minh City, Vietnam [KHV] - HG02064, HG02067

  • ○ Japanese in Tokyo, Japan (HapMap) [JPT] - NA19068*

• Africa (AFR):

  • ○ Luhya in Webuye, Kenya (HapMap) [LWK] - GM19440*

  • ○ Gambian in Western Division, The Gambia [GWD] - HG02623*

  • ○ Esan from Nigeria [ESN] - HG03115*

• Europe (EUR):

  • ○Toscani in Italia (Tuscans in Italy) (HapMap) [TSI] - GM20587*

  • ○ British from England and Scotland, UK [GBR] - HG00250*

  • ○ Finnish in Finland [FIN] - HG00353*

• America (AMR):

  • ○ Mexican Ancestry in Los Angeles, California, USA (HapMap) [MXL] - GM19789*

  • ○ Peruvian in Lima, Peru [PEL] - HG01971*

Asterisks next to sample IDs refer to samples processed at UCSF. ~ 99% of the SNPs are phased in all the samples. For all the sequences, < 1% of each sequence has zero coverage. There are, however, differences in the exact protocols used for the samples sequenced at Genentech and UCSF. As a result, lengths of the phase blocks as well as the N50 values for the phase blocks differ by a factor of 10 between the two sets of samples. However, even the smallest phase blocks are long enough for accurate phasing. Statistics for the experimental sequencing like sequence coverage, N50, and fraction of SNPs phased can be found in the Additional file 2.

Preprocessing 1000 genomes data

The 1000 Genomes data was separated into individual and chromosome specific VCFs using vcftools [27]. Further, the variants were filtered for biallelic SNPs, phased (i.e. variants already phased in the 1000 Genomes VCFs [8]), filtered for PASS, and indels were removed. The experimentally phased data also had a very small fraction of unphased SNPs, which were removed by filtering with vcftools. The analysis was performed only for autosomes.

Phasing analysis

The alternate (ALT) allele frequencies of all the SNPs of interest were obtained from the 1000 Genomes data and converted to minor allele frequencies to be able to analyze switch error as a function of minor allele frequencies. The filtered SNPs from the experimental data were split into phase sets, based on phase set information available in the experimental VCF files. Long runs of homozygosity, leading to uncertainty in the phasing method associated with the experimental sequencing cause the phasing to be broken off. This leads to the creation of multiple phase sets in the final experimental sequences [28]. Switch error was calculated between the experimental and 1000 Genomes data for each phase set in each chromosome of each individual from the experimental dataset. Switch error is defined as percentage of possible switches in haplotype orientation used to recover the correct phase in an individual [29] or equivalently, proportion of heterozygous positions whose phase is wrongly inferred relative to the previous heterozygous position [30]. vcftools returns the switch error as well as all positions of switches occurring along the chromosome.

Switch error as a function of minor allele frequency

ALT allele frequencies were accessed for each of the switch positions (i.e. both heterozygotes at the ends of each out-of-phase segment) from the data and were converted to minor allele frequencies. A distribution of all the switch positions as a function of minor allele frequency was plotted for each chromosome in each individual.

Switch error as a function of inter SNP distance

Positions of each SNP were accessed from the data. The number of intermediate switches were counted for all pair of SNPs, not only consecutive SNPs. If the number of switches between two SNPs were odd, a switch error was counted. This was used to calculate the distribution of switch errors as a function of inter-SNP distance.

Imputation analysis

The entire imputation analysis is performed for each chromosome for each individual.

Generate recombination map

IMPUTE v2 [11] makes available recombination maps for each chromosome using the 1000 Genomes data for the GRCh37 assembly. A recombination map was obtained for each chromosome for GRCh38 by lifting over the GRCh37 maps using the liftOver [31] software. ~ 8 k positions (0.2%) were removed from the lifted over recombination map because liftover resulted in them being in the incorrect order.

Generate reference panel

A reference haplotype panel was generated for all individuals from the 1000 Genomes data by subsetting it to the specific population of interest. 1000 Genomes data for the individuals which were experimentally sequenced was not included in the reference panel. Vcftools was used to filter out the individuals of interest from the 1000 Genomes data. Bcftools was used to convert the VCF data to haps-sample-legend format. An alternate approach was also used, where the entire 1000 Genomes data was used to generate a reference haplotype panel. The number of haplotypes in the population specific reference panels were: AFR-1316, AMR-690, EUR-1000, EAS-990, SAS-956.

Generate study panel

A study panel was generated for the experimentally sequenced individuals selected. The study panel is assumed to be genotyped at positions corresponding to the Illumina InfiniumOmni2.5–8 array. Array positions were lifted over from GRCh37 to GRCh38 using liftOver. 1000 Genomes haplotypes (since 1000 Genomes data is prephased, the study panel is also in the form of haplotypes rather than genotypes) from the 1000 Genomes final calls for those positions for those individuals were selected to create the study panel using vcftools. Filtered VCF files were converted to the haps-sample format using bcftools.

Run imputation

Missing positions are imputed using IMPUTE v2. Imputation was performed in 5 Mb windows. The genotype output by imputation was converted to VCF format using bcftools. VCFs produced over all windows were combined using vcf-concat. IMPUTE v2 generally phases the typed genotyped sites in study panel. This is followed by imputation. IMPUTE v2 then performs an iterative process performing multiple Monte-Carlo steps alternating phasing and imputation. For this analysis, however, as haplotypes from the 1000 Genomes project were directly used to generate the study panel, the phasing step was not performed.

Filter positions

For one part of the analysis, i.e. estimating errors in the positions represented in the experimentally phased VCFs (called experimental SNPs throughout the manuscript), the positions from those VCFs were filtered from the imputed data using vcftools. Experimental genotypes from the experimental VCFs were obtained for each individual of interest using vcftools. SNPs with duplicate entries in either the imputed or experimental data were removed. Continent-specific allele frequencies were obtained for the experimental SNPs from the 1000 Genomes data using vcftools, to be able to analyze switch error as a function of Minor Allele Frequencies. For the other part of the analysis, i.e. estimating errors for all positions in the 1000 Genomes data, the allele fractions were similarly obtained for all of the SNPs.

Imputation error

Imputation error was computed in the form of genotype discordance (fraction of genotypes being incorrectly identified). Imputation error was computed for both, the SNPs in the experimental data and all the SNPs in 1000 Genomes data. Error is computed as a function of minor allele frequency. The continent-specific minor allele frequencies were used for analyzing the imputation error. r² between the imputed and experimental genotypes for each SNP is another common method used to estimate imputation accuracy, and is considered to minimize the dependence on the allele frequency. However, we only have between 2 and 11 individuals in each continental group experimentally sequenced and phased in our experiments, which are too small numbers to be able to compute an r² value for each SNP. Hence r² values have been computed for all SNPs in each allele frequency window. These windows are computed with alternative allele frequencies instead of minor allele frequencies to allow comparison with previously estimated imputation accuracies [3].

For all analysis where error rate is computed as a function of the continent-specific minor allele frequency (genotyping error and imputation error; Figs. 1, 2, 7, 8), the minor allele frequencies are binned as MAF = 0.0%, 0.0–0.2%, 0.2–0.5%, 0.5–1%, 1–5%, MAF > = 5%. For the analysis where all 1000 Genomes minor allele frequencies are used (phasing error and imputation error comparing use of multiple reference panels; Figs. 3, 4, 10), the minor allele frequencies are binned into only five bins, i.e. there is no MAF = 0.0% bin. Rest of the bins are the same as for the continent-specific MAF bins.

In addition, imputation accuracy was also computed as genotype r² (Fig. 9) for all 1000GP SNPs. This is plotted against alternate allele frequency (instead of minor allele frequency) to enable comparison with the previous accuracy estimates in the 1000GP phase 3 paper [3]. r² values are computed for all genotypes values of all SNPs in each alternative allele frequency (AAF) bin instead of per SNP to deal with the fact that the AFR, AMR, and EUR populations have only 3, 2, and 3 individuals respectively. The AAF values are binned as AAF < 0.2%, 0.2–0.5%, 0.5–1%, 1–2%, 2–5%, 5–10%, 10–20%, 20–50% and 50–100%.

Experimental methods

Samples processing

High Molecular Weight (HMW) Genomic DNA was extracted and converted into 10x sequencing libraries according to the 10X Genomics (Pleasanton, CA, USA) Chromium Genome User Guide and as published previously [28]. Briefly, Gel Bead-in-Emulsions (GEMs) were made with 1.25 ng HMW template gDNA, Master-mix Genome Gel Beads and partitioning oil on the microfluidic Genome Chip. Isothermal incubation of the GEMs (for 3 h at 30 °C; for 10 min at 65 °C; stored at 4 °C) produced barcoded fragments ranging from a few to several hundred base pairs. After dissolution of the Genome Gel Bead in the GEM Illumina Read 1 sequencing primer, 16 bp 10x barcode and 6 bp random primer are released. The GEMs were then broken and the pooled fractions were recovered. Silane and Solid Phase Reversible Immobilization (SPRI) beads were used to purify and size select the fragments for library preparation. Library prep was performed according to the manufacturer’s instructions described in the Chromium Genome User Guide Rev. C. Libraries were made using 10x Genomics adapters. The final libraries contain the P5 and P7 primers used in Illumina bridge amplification. The barcoded libraries were then quantified by qPCR (KAPA Biosystems Library Quantification Kit for Illumina platforms). Sequencing was done using Illumina HiSeq 4000 with 2 × 150 paired-end reads. Raw reads were processed, aligned to the reference genome, and had SNPs called and phased using 10X Genomics’ Long Ranger software (version 2.1.1 or 2.1.6) with the “wgs” pipeline with default settings.

The sequence data for individuals sequenced experimentally can be found in BioProjects PRJNA418343 and PRJNA435626 (samples sequenced at UCSF) and BioProject PRJNA544309 (samples sequenced at Genentech).

1000GP:

1000 Genomes Project

AAF:

Alternate Allele Frequency

AFR:

Africa

ALT:

Alternate [Alleles]

AMR:

America

BEB:

Bengali in Bangladesh

CDX:

Chinese Dai in Xishuangbanna, China

CHB:

Han Chinese in Beijing, China

EAS:

East Asia

ESN:

Esan from Nigeria

EUR:

Europe

FIN:

Finnish in Finland

GBR:

British from England and Scotland

GEMs:

Gel Bead-in-Emulsions

GIH:

Gujarati Indians in Houston, Texas

GWAS:

Genome Wide Association Study

GWD:

Gambian in Western Division, The Gambia

HMW:

High Molecular Weight

HRC:

Haplotype Reference Consortium

IBD:

Identical-By-Descent

ITU:

Indian Telugu in the UK

JPT:

Japanese in Tokyo, Japan

KHV:

Kinh in Ho Chi Minh City, Vietnam

LWK:

Luhya in Webuye, Kenya

MAF:

Minor Allele Frequency

MXL:

Mexican Ancestry in Los Angeles, California, USA

PEL:

Peruvian in Lima, Peru

PJL:

Punjabi in Lahore, Pakistan

qPCR:

quantitative Polymerase Chain Reaction

SAS:

South Asia

SNP:

Single Nucleotide Polymorphism

SPRI:

Silane and Solid Phase Reversible Immobilization

STU:

Sri Lankan Tamil in the UK

TSI:

Toscani in Italia

VCF:

Variant Call Format

WGS:

Whole Genome Sequencing

Not applicable.

Funding

JDW was supported in part by NIH grant R01 GM115433. SB was supported by Genentech research grant CA0095684. Neither funding body had any influence over study design, analysis and interpretation of the data, or in writing the manuscript.

Affiliations

1. Institute for Human Genetics, University of California, San Francisco, CA, 94143, USA

   • Saurabh Belsare
   • , Pui-Yan Kwok
   •  & Jeffrey D. Wall

2. Department of Dermatology, University of California, San Francisco, CA, 94143, USA

   • Michal Levy-Sakin
   • , Yulia Mostovoy
   •  & Pui-Yan Kwok

3. Department of Molecular Biology, Genentech Inc, 1 DNA Way, South San Francisco, CA, 94080, USA

   • Steffen Durinck
   • , Subhra Chaudhuri
   • , Andrew S. Peterson
   •  & Somasekar Seshagiri

4. School of Biomedical Science, Engineering, and Health Systems, Drexel University, Philadelphia, PA, 19104, USA

   • Ming Xiao

5. Cardiovascular Research Institute, San Francisco, CA, 94143, USA

   • Pui-Yan Kwok

6. Department of Epidemiology and Biostatistics, University of California, San Francisco, CA, 94143, USA

   • Jeffrey D. Wall

Contributions

SB contributed to the design of and performed the computational analysis and was a major contributor in writing the manuscript. MSL, YM, SD, SC, MX contributed to the design of and performed the experimental sequencing of the individual sequences. ASP, PYK, and SS contributed to the design of the experimental sequencing of the individual sequences. JDW contributed to the design of the computational analysis and was a major contributor in writing the manuscript. All authors have read and approved the final manuscript.

Corresponding authors

Correspondence to Saurabh Belsare or Jeffrey D. Wall.

Ethics approval

Not applicable.

Consent for publication

Not applicable.

Competing interests

Genentech authors hold shares in Roche. The other authors declare no conflicts of interest.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Reprints and Permissions