Source-linked AI summary
Integrating sequencing datasets to form highly confident SNP and indel genotype calls for a whole human genome
Justin M. Zook, Brad Chapman, Jason Wang, David Mittelman, Oliver Hofmann, Winston Hide, Marc Salit
TL;DR
Accurate whole-genome genotype benchmarks are needed to assess sequencing accuracy, but methods and algorithms can produce discordant calls. This paper integrates multiple sequencing datasets to generate highly confident NA12878 calls that are more sensitive, specific, and less biased than individual datasets and publicly available for benchmarking.
Problem
Discordant sequencing methods and algorithms created a need for highly accurate whole-genome genotype calls for benchmarking.
Method
The authors integrate and arbitrate genotype evidence across multiple sequencing technologies, datasets, mappers, and callers to identify confident and uncertain calls.
Results
The integrated calls are more sensitive, specific, and less biased than individual datasets, with 564,410 SNPs correctly called on a microarray.
Takeaways & Limitations
The publicly available calls provide a whole-genome resource for assessing biases and genotype-call accuracy across sequencing and bioinformatics methods.
Takeaways & Limitations
Comparisons to the integrated genotypes exclude approximately 23% of the genome, including potential structural variants and other difficult regions.
Abstract
from arXiv · showhide
Clinical adoption of human genome sequencing requires methods with known accuracy of genotype calls at millions or billions of positions across a genome. Previous work showing discordance amongst sequencing methods and algorithms has made clear the need for a highly accurate set of genotypes across a whole genome that could be used as a benchmark. We present methods to make highly confident SNP, indel, and homozygous reference genotype calls for NA12878, the pilot genome for the Genome in a Bottle Consortium. We minimize bias towards any method by integrating and arbitrating between 14 datasets from 5 sequencing technologies, 7 mappers, and 3 variant callers. Regions for which no confident genotype call could be made are identified as uncertain, and classified into different reasons for uncertainty. Our highly confident genotype calls are publicly available on the Genome Comparison and Analytic Testing (GCAT) website to enable real-time benchmarking of any method.
Different variant representations make comparison difficult
Indels and complex variants can be represented correctly in multiple ways, making comparisons across variant callers difficult. The study regularized variant call files to avoid counting equivalent representations as different variants.
- Indels and complex variants are difficult to compare across callers because the same variant can have multiple correct representations.Nearby SNPs and indels may be expressed in different but equivalent forms.
- The variant call files were regularized with vcflib’s vcfallelicprimitives module to prevent equivalent variant expressions from being counted separately.Regularization minimized counting nearby SNPs and indels expressed differently as distinct variants.
- A complex variant was represented four different ways by datasets using different mappers, despite replacing the same 6 bases with 5 bases.The example involved reads from Illumina, 454, and Complete Genomics datasets mapped with BWA, ssaha2, CGTools, and Novoalign.
Integrated variant calls are highly sensitive and specific
Integrated genotype calls showed high sensitivity and specificity across orthogonal microarray and sequencing benchmarks, while their elevated Ti/Tv for novel whole-genome variants indicated a lower error rate. Comparisons also identified limitations of microarrays and highlighted uncertain genomic regions requiring further inspection.
- Call accuracy: Integrated calls had a higher Ti/Tv than other datasets for novel whole-genome variants, which usually indicates a lower error rate.For common variants, integrated-call Ti/Tv was comparable to other datasets in both whole genome and exome analyses.
- Orthogonal validation: 564,410 microarray SNPs were correctly called, while 1,332 microarray SNPs were absent and 527 variants in the integrated set were at homozygous-reference microarray positions.The comparison assessed sequencing calls at sites interrogated by the microarray.
- Orthogonal validation: 427 SNPs and 42 indels from GeT-RM were correctly genotyped, and calls were concordant for all 124 SNPs and 37 indels in XPrize Sanger data.These comparisons used high-quality variants found across multiple sequencing platforms and independent Sanger sequencing data.
- Call accuracy: 208 freebayes variants with coverage greater than 20 were called highly confident homozygous reference, and all manually inspected random 10 % appeared likely false positives.This comparison was used to test whether the integrated methods missed sites identified by another variant caller.
- Benchmark limitations: Microarray-based sensitivity assessment is limited because microarrays cover only uniquely probe-accessible regions, whereas integrated calls contain a 67 times higher percentage of low-complexity regions.The passage notes that 9.8% of uncertain sites have low complexity, motivating whole-genome benchmarking beyond microarray-accessible regions.
- Uncertain regions: About 1 million variants in the 250 bp Illumina HaplotypeCaller VCF lie inside NA12878 structural variants reported to dbVar, with over 200,000 variants in several other uncertain categories.Those categories include unresolved conflicting genotypes, segmental duplications, low coverage or mapping quality, and simple repeats.
Discussion
The study establishes methods for integrating multiple sequencing datasets into highly confident whole-genome SNP and indel genotype calls that reduce method-specific bias. It also emphasizes that benchmark construction and interpretation affect performance assessment, motivating careful review and ongoing refinement of the call set.
- Contributions: Integrating datasets across sequencing technologies produces genotype calls more sensitive, specific, and less biased than any individual dataset.The methods account for systematic sequencing, local alignment, and mapping errors associated with individual datasets.
- Benchmark limitations: Excluding difficult genomic regions or labeling them uncertain can lower the apparent false-positive and false-negative rates of assessed datasets.The passage warns that benchmark construction changes performance assessment when difficult regions are down-weighted or disregarded.
- Benchmark limitations: Benchmark calls should have much lower FP and FN rates than the corresponding error rates of the assessed dataset.The authors compare integrated calls with independent methods to evaluate bias and error rates.
- Use and refinement: Users should inspect alignments around discordant genotype calls because overall sensitivity and specificity may not resolve questionable variants.Manual inspection across multiple datasets can help identify difficult regions and questionable variants.
- Clinical translation: The methods provide a basis for reference materials, data, and methods intended to support translating genome sequencing into clinical practice.The Genome in a Bottle Consortium will develop these resources, while the call set will be maintained and refined as new methods become available.
Online Methods
Integrated genotype calls for Coriell DNA sample NA12878 used nine whole-genome and three exome sequencing datasets. These datasets spanned five sequencing platforms, including six Illumina whole-genome and two Illumina exome datasets.
- Dataset integration: 9 whole-genome and 3 exome sequencing datasets were integrated to form genotype calls for Coriell DNA sample NA12878.Table 1 describes each dataset’s source, platform, mapping algorithm, coverage, and aligned read length.
- Sequencing platforms: 6 whole-genome and 2 exome datasets came from Illumina sequencers, while SOLiD, 454, Complete Genomics, and Ion Torrent contributed additional datasets.The additional contributions comprised one whole genome each from SOLiD, 454, and Complete Genomics, plus one exome from Ion Torrent.
- Dataset composition: The whole-genome datasets included 2 PCR-free datasets among the 6 generated on Illumina sequencers.Some datasets had BAM files published.
Comparison of variant calls using different methods
Variant comparisons first normalized short indels, complex variants, and MNPs so equivalent representations across methods were not counted as discordant. The scripts and resulting BED file were made publicly available.
- Variant representation was regularized for short indels, complex variants, and MNPs to prevent equivalent method-specific expressions from being counted as discordant.The workflow used vcflib’s vcfallelicprimitives module to help normalize representations.
- The Perl scripts and resulting BED file were publicly released on GitHub.
Obtaining highly confident genotypes for training VQSR
Candidate SNP and indel sites were gathered across datasets, then genotypes were forcibly recalled and combined to produce highly concordant calls. These calls independently trained VQSR models by genotype class using predominantly single-Gaussian annotation distributions.
- Candidate site discovery: A union of low-threshold GATK calls, Complete Genomics SNP calls, and Cortex assembly calls identified possible SNP and indel sites across datasets.GATK UnifiedGenotyper and HaplotypeCaller used a variant quality score threshold of 2; Complete Genomics used its unfiltered CGTools 2.0 SNP set.
- Genotype recall: Genotypes were forcibly recalled for every dataset at candidate SNPs, while HaplotypeCaller performed local de novo assembly around candidate indels and complex variants.UnifiedGenotyper recalled candidate SNP genotypes individually, and HaplotypeCaller recalled indels and complex variants individually before combining calls.
- Genotype recall: Combined UnifiedGenotyper and HaplotypeCaller calls prioritized HaplotypeCaller within 20bp of a HaplotypeCaller indel.The recalled calls were combined with preference given to HaplotypeCaller in the specified indel-proximal region.
- VQSR training: Highly concordant heterozygous and homozygous variant genotypes independently trained VQSR Gaussian Mixture Models for each dataset and genotype class.Training was separated for heterozygous and homozygous variant and reference genotypes because their annotation distributions and bias characteristics could differ.
- VQSR training: Each annotation was modeled with a single Gaussian because additional Gaussians often fit noise and frequently prevented model convergence.The single-Gaussian choice reflected approximately Gaussian annotation distributions and was applied to annotations associated with alignment, mapping, and systematic sequencing errors.
Arbitration between datasets with conflicting genotypes
Conflicting genotype datasets are arbitrated by favoring calls from datasets with typical site-specific annotations and requiring concordance across multiple datasets. The method also flags difficult regions and applies tranche-based quality requirements to reduce incorrect calls.
- Arbitration between datasets with conflicting genotypes: Datasets with unusual site-specific annotations associated with bias are considered less likely to be correct than datasets with typical annotations.For each possible variant site, the approach first determines whether at least two datasets confidently call the same genotype.
- Arbitration between datasets with conflicting genotypes: Regions difficult for all methods remain uncertain even when datasets agree, including cases involving reference-assembly copy-number differences and homozygous SNPs.VQSR tranches identify annotations associated with systematic sequencing errors, alignment problems, mapping problems, and atypical allele balance.
- Arbitration between datasets with conflicting genotypes: For homozygous reference genotypes, at least 2 datasets must have an alignment tranche < 99; heterozygous genotypes require multiple tranche thresholds below 95 or 99.Heterozygous calls require at least 3 datasets with mapping tranche < 99, plus at least 2 datasets each meeting systematic sequencing error, alignment, mapping, and allele balance tranche < 95.
GCAT performance assessment of dataset
GCAT comparisons used regularized, structurally filtered variant files to assess integrated calls against fosmid calls. Manual curation resolved discordances and suggested very low error rates for the integrated calls in highly confident regions.
- Comparison preparation: GCAT comparisons first regularized variants with vcflib vcfallelicprimitives and excluded dbVar structural variants from whole-genome comparisons.Whole-exome comparisons additionally used the manufacturer’s target exome bed file.
- Discordance resolution: Almost all variants present in the integrated calls but absent from fosmid calls were fosmid false negatives caused by mis-called complex variants or overly stringent filtering.One SNP, chr2:108078636, was clearly homozygous variant in all integrated datasets but homozygous reference in the fosmids, with the reason unclear.
- Performance estimate: ~3 partial complex variant calls and between 0 and 1 false positive or false negative simple SNP or indel calls per 30 million highly confident bases were estimated for the integrated calls.The estimate was based on 119 million fosmid reference bases, 25% manual curation, ~94,500 TP SNPs, and ~1400 TP indels.