Accurate and comprehensive sequencing of personal genomes

Abstract

As whole-genome sequencing becomes commoditized and we begin to sequence and analyze personal genomes for clinical and diagnostic purposes, it is necessary to understand what constitutes a complete sequencing experiment for determining genotypes and detecting single-nucleotide variants. Here, we show that the current recommendation of ∼30× coverage is not adequate to produce genotype calls across a large fraction of the genome with acceptably low error rates. Our results are based on analyses of a clinical sample sequenced on two related Illumina platforms, GAII(x) and HiSeq 2000, to a very high depth (126×). We used these data to establish genotype-calling filters that dramatically increase accuracy. We also empirically determined how the callable portion of the genome varies as a function of the amount of sequence data used. These results help provide a "sequencing guide" for future whole-genome sequencing decisions and metrics by which coverage statistics should be reported.

Figure 1.

Breadth versus depth of whole-genome coverage. The x-axis represents the minimum number of high-quality bases (≥Q20) from high-quality alignments (≥MapQ30), and the y-axis represents the proportion of genome (A) or coding exome (B) covered at that depth. To calculate percentages, the total size of hg18 build and the total number of non-redundant coding bases from the UCSC Known Genes table (2,852,680,119 bp and 34,068,542 bp, respectively) were used. Gaps and pseudo-autosomal regions (PAR) were excluded. Values were plotted for GAII_x (triangle), HiSeq flowcell A (orange square), HiSeq 2000 flowcell B (dark red square), and all data sets combined (circle).

Figure 2.

Effect of alignment filter on the discordance rate of identical genomes. The number of discordant positions (y-axis) was observed by varying MapQ values (x-axis). A MapQ value of 0 indicates that no mapping quality filter was applied.

Figure 3.

Determination of genotype confidence threshold for genotype calls. (A,B) The x-axes represent Q20 depth for genotype calls from one of the 50× genomes, and the y-axes represent corresponding MPG scores. (A) A random set of ∼8700 concordant genotypes; (B) 8710 discordant genotypes. Black lines represent a line with slope of 0.5, which is the confidence threshold used to filter genotypes. (C) The fraction of genotypes retained by varying the confidence threshold; (blue curve) the fraction of concordant genotypes retained; (red curve) the fraction of discordant genotypes retained.

Figure 4.

Comparison of identical genomes at various mapped depths. (A,B) The x-axes represent the average mapped depths at which two identical genomes were compared. The y-axes represent the proportion of hg18 callable in both genomes (A) and the discordance per megabase of callable sequence (B). Analyses were done on all unique alignments (MapQ > 0) without applying any filters (red curve) and after applying mapping quality and genotype confidence filters as explained in the text (blue curve).

Figure 5.

Genotype calling as a function of average mapped depth. The x-axes represent the average mapped depth of each data set, and the y-axes represent the proportion of the whole genome (dark blue circles) and coding exome (green triangles) that is callable (A), the number of SNVs detected (B), the proportion of Illumina BeadChip positions callable (C), and the concordance rates with the BeadChip calls (D).

Figure 6.

Improved representation of genome with TruSeq v3 sequencing chemistry and software. The x-axis represents the average mapped depth of the data set, and the y-axis represents the proportion of the whole genome (dark blue circles) and coding exome (green triangles) that is callable.