PaperlayerView exact revision on arXiv

Source-linked AI summary

Towards Better Understanding of Artifacts in Variant Calling from High-Coverage Samples

Heng Li

arXiv:1404.0929v2q-bio.GN

TL;DR

Variant-calling error rates and causal artifacts remain difficult to measure without an unbiased whole-genome truth set. The paper uses haploid CHM1 and diploid NA12878 data across multiple mappers and callers to investigate false heterozygous calls. It identifies low-complexity-region realignment and incomplete reference representation as major error sources, estimating raw errors at 1 in 10–15kb and post-filtered errors at 1 in 100–200kb without significant sensitivity loss.

  • Problem

    The whole-genome error rate of SNP and INDEL calls and the leading artifacts remained unclear because existing evaluations lacked an unbiased truth set and often hid systematic errors.

  • Method

    The study compares multiple mapper–caller call sets from haploid CHM1 and diploid NA12878, using CHM1 heterozygous calls to investigate errors and test filters.

  • Results

    1 in 10–15kb raw genotype calls are errors, while post-filtered call sets reach approximately 1 in 100–200kb without significant compromise on sensitivity.

  • Takeaways & Limitations

    Erroneous realignment in low-complexity regions and incomplete reference representation are major error sources, supporting continued improvements in both areas.

  • Takeaways & Limitations

    The maximum depth filter applies to high-coverage data with uniform read depth and is ineffective for exome or shallow-coverage data.

Abstract

from arXiv · show

Motivation: Whole-genome high-coverage sequencing has been widely used for personal and cancer genomics as well as in various research areas. However, in the lack of an unbiased whole-genome truth set, the global error rate of variant calls and the leading causal artifacts still remain unclear even given the great efforts in the evaluation of variant calling methods. Results: We made ten SNP and INDEL call sets with two read mappers and five variant callers, both on a haploid human genome and a diploid genome at a similar coverage. By investigating false heterozygous calls in the haploid genome, we identified the erroneous realignment in low-complexity regions and the incomplete reference genome with respect to the sample as the two major sources of errors, which press for continued improvements in these two areas. We estimated that the error rate of raw genotype calls is as high as 1 in 10-15kb, but the error rate of post-filtered calls is reduced to 1 in 100-200kb without significant compromise on the sensitivity. Availability: BWA-MEM alignment: http://bit.ly/1g8XqRt; Scripts: https://github.com/lh3/varcmp; Additional data: https://figshare.com/articles/Towards_better_understanding_of_artifacts_in_variating_calling_from_high_coverage_samples/981073

1 INTRODUCTION

The paper addresses the unresolved whole-genome error rate of SNP and INDEL calling and the difficulty of identifying systematic artifacts without an unbiased truth set. It evaluates variant-calling accuracy using a haploid human cell line, where most heterozygous calls are expected to be errors.

  • Limits of existing evaluation: Simulation-based evaluations cannot reproduce artifacts from non-random variants, dependent errors, incomplete references, copy-number variation, or large-scale effects.Genome-to-genome simulation is also limited for mammalian genomes because references are incomplete and whole-genome alignment is imperfect.
  • Limits of existing evaluation: Experimental validation and array-based comparisons are biased toward easier genomic regions and may themselves contain artifacts or errors.Comparing call sets measures relative rather than absolute accuracy when pipelines share artifacts.
  • Limits of existing evaluation: Pedigree and same-individual comparisons can identify many inconsistencies as errors but are insensitive to errors recurring across all samples.Systematic errors at a locus remain unidentified when every sample receives the same erroneous heterozygous call.
  • Research gap: The global error rate and leading sources of SNP and INDEL calling errors remained unclear, especially for systematic errors.Existing evaluations often hid systematic errors or relied on estimates from easy sites.
  • Study approach: The study uses a haploid human cell line to evaluate variant calling, treating most heterozygous calls as errors across the genome.It combines multiple read mappers and variant callers with caller-oblivious genotyping and filtering.

2 DATA SETS AND DATA ANALYSIS

The study analyzes deep Illumina data from haploid CHM1 and diploid NA12878 cell lines using multiple alignments and variant callers. It standardizes representation and genotype evaluation while documenting coverage, duplicate handling, and preprocessing choices.

  • Data sets: CHM1 is haploid, so heterozygous variant calls are treated as errors, while NA12878 provides a diploid control for sensitivity.The study focuses on deep sequencing data from two cell lines.
  • Data sets: CHM1 uses six runs from one library at about 65-fold coverage before duplicate removal, whereas NA12878 paired-end data provide about 55-fold coverage without PCR amplification.The full CHM1 data set exceeds 100-fold coverage, but only six runs are analyzed.
  • Alignment: Reads are mapped with Bowtie2 and BWA-MEM for CHM1, and with BWA, Bowtie2, and BWA-MEM for NA12878.Most analyses map reads to the hs37d5 reference genome.
  • Preprocessing: Picard marks 20% of CHM1 reads as PCR duplicates, while its 1.5% estimate for NA12878 is treated as false positive and excluded from subsequent analysis.NA12878’s library was constructed without amplification.
  • Variant calling: The evaluation calls SNPs and short INDELs with FreeBayes, GATK UnifiedGenotyper, Platypus, SAMtools, and GATK HaplotypeCaller.Complex variants are decomposed into individual events, and overlapping variants are reduced to the higher-quality event where applicable.
  • Genotyping: Genotypes are re-called from caller-provided genotype likelihoods to reduce differences caused by callers’ distinct genotype priors.Platypus-reported genotypes are retained for multiallelic variants because it does not provide genotype likelihoods for them.

2.4 Variant filtering

The study applies universal filters across call sets and estimates accuracy by comparing heterozygous calls in haploid CHM1 with calls in diploid NA12878. The filters target genomic complexity, depth, allele balance, strand support, and quality.

  • Filtering strategy: Universal filters are chosen to avoid caller-specific filtering that would complicate comparisons and obscure artifacts.The filter set includes low complexity, maximum depth, allele balance, double strand, Fisher strand, and quality filters.
  • Filtering strategy: The low-complexity filter removes variants overlapping mdust-identified low-complexity regions, which cover 2.0% of autosomal bases in GRCh37.This regional filter is independent of alignment and variant calling.
  • Filtering strategy: The maximum depth filter removes sites covered by an excessive number of reads, while allele balance removes sites with too few non-reference reads.Depth definitions can differ between callers.
  • Filtering strategy: Double-strand and Fisher-strand filters target strand support and reference/non-reference strand correlation, respectively.The double-strand filter is unavailable for GATK calls because GATK does not report the required strand counts.
  • Accuracy measurement: The CHM1-to-NA12878 comparison estimates heterozygous false-positive rate as Nh/Nd and sensitivity as Nd − Nh.Nh is the number of heterozygous CHM1 calls and Nd is the number of NA12878 calls.
  • Accuracy measurement: The study manually reviews more than 200 heterozygous INDELs using alignment visualization, flanking-region inspection, and read assembly to characterize error modes.The review examines alignment quality, positional biases, and reference-genome complexity.

3 RESULTS

The results present filters sequentially to clarify their effects across call sets. Under a possible somatic-mutation interpretation of CHM1 heterozygotes, call-set error rates are approximately 1 in 100–200kb.

  • Filters were initially applied independently to each call set, then presented in an ordered sequence for clarity.Later filters depended on those applied earlier in the presented analysis.
  • Figure 1 overviews the breakdown of various filters across multiple call sets.
  • 1 in 100–200kb: call sets generally have this approximate error rate when CHM1 may contain true heterozygotes from somatic mutations.The passage equates this rate with approximately 15,000–30,000 false heterozygous calls, but the sentence is truncated.

3.1 Checking the ploidy of CHM1

The small overlap among CHM1 heterozygous call sets supports the conclusion that most heterozygotes arise from sequencing or algorithmic artifacts rather than non-haploid biology. Even conservatively, the estimated error rate reaches one per 170kbp.

  • Establishing CHM1 ploidy: Most CHM1 heterozygotes are likely sequencing or calling artifacts because different callers identify largely different subsets of errors.A genuinely non-haploid sample would produce many heterozygotes shared across callers, whereas artifacts produce low-consistency call sets.
  • Establishing CHM1 ploidy: Manual review found that more than half of heterozygous SNPs called by multiple callers appear real across CHM1 and additional Illumina libraries.These SNPs generally had strong alignment support, average read depth, and no overlap with known segmental duplications or error-prone motifs.
  • Filtering analysis: The analysis applied filters sequentially, with maximum depth after low-complexity filtering and miscellaneous filters after both earlier stages.The miscellaneous filters used allele balance, strand support, and Fisher strand criteria.
  • Estimating the error rate: 1 error per 170kbp is an upper-bound estimate even when all intersecting calls are assumed to be genuine CHM1 variants.The estimate is calculated as 3Gbp divided by 17.7k.
  • Estimating the error rate: 10 times more raw heterozygous calls occur in NA12878 than in CHM1, supporting the interpretation that CHM1 heterozygotes largely reflect major sequencing or calling artifacts.This comparison is used alongside the conservative CHM1 error-rate estimate.

3.2 The low-complexity filter

Low-complexity regions occupy only 2% of the human genome but contain most CHM1 heterozygous INDEL calls and many heterozygous SNPs. Manual review links these errors to PCR effects and faulty realignment, while local assembly can resolve difficult cases.

  • Low-complexity enrichment: 2% of the human genome contains 80–90% of CHM1 heterozygous INDEL calls and up to 60% of heterozygous SNPs.These low-complexity regions therefore concentrate a disproportionate share of apparent heterozygous calls.
  • Low-complexity enrichment: The false-positive rate in low-complexity regions ranges from 10% to 40%, compared with about 1–8% for heterozygous INDELs outside them.The ranges vary depending on the call set.
  • Error mechanisms: Two major sources of INDEL genotyping errors in low-complexity regions are potential PCR errors and realignment errors.Manual review of over 100 sites identified both mechanisms.
  • Error mechanisms: PCR-free NA12878 data still show smaller call-set intersections in low-complexity than high-complexity regions, indicating additional non-PCR error sources.PCR errors are believed to affect SNPs to a lesser extent, so the remaining disagreement points to other mechanisms.
  • Realignment errors: About 50–70% of reviewed heterozygous INDELs larger than 1bp can be corrected with better realignment.Variant callers often fail to produce the optimal realignment in low-complexity regions.
  • Realignment errors: Local assembly with fermi is frequently more effective than INDEL callers and can make difficult realignment errors trivial to resolve.The proposed advantages include reference independence, long-range consistency, and topology-based error cleaning.

3.3 The maximum read depth filter

The maximum-depth filter is the most effective evaluated filter against false heterozygotes, particularly for BWA-MEM alignments, while alignment and reference-related artifacts explain important residual errors.

  • Figure 5 uses CHM1 heterozygous SNP counts as a false-positive proxy and the CHM1–NA12878 count difference as a sensitivity proxy.Curves nearer the top-left indicate better error classification.
  • The maximum-depth filter is most effective against false heterozygotes, especially those from BWA-MEM alignments.At mean depth d ≈50, thresholds between d + 3d and d + 4d remove many false positives with little sensitivity loss.
  • Many false positives removed by maximum-depth filtering are attributed to CNVs or paralogous sequences absent from the human reference genome.These artifacts create excessive read depth and mismapping signals.
  • Local misalignments can generate multiple heterozygotes around low-complexity INDELs, whereas local assembly can resolve difficult examples.Bowtie2 aligned none of the example reads with gaps, BWA-MEM aligned four correctly, and HaplotypeCaller locally assembled reads.
  • Conservative Bowtie2 mapping-quality estimates reduce susceptibility to CNV and incomplete-reference artifacts but may increase false negatives.One example produced a SNP from BWA-MEM but not Bowtie2; long-INDEL call differences between mappers remained unresolved.
  • The maximum-depth filter applies mainly to high-coverage data with uniform read depth and is ineffective for exome or shallow-coverage sequencing.A Hardy–Weinberg filter was nearly as effective on CHM1 but required multiple population-, mapper-, and caller-dependent thresholds.

3.4 Other filters

Additional filters remove some false heterozygotes, but their effectiveness and optimal thresholds vary across call sets and trade off false-positive against false-negative rates.

  • AB, DS, FS, and QU filters remove additional CHM1 false heterozygotes, but their effectiveness varies across call sets.The authors selected reasonable thresholds from ROC-like curves rather than optimizing separately for every call set.

3.5 Effect of PCR duplicates

PCR duplicates affected unfiltered heterozygous calls, especially INDELs, but filtering produced similar total SNP and INDEL counts with or without duplicates.

  • 20% of CHM1 data were discarded because of PCR duplicates.Calling without duplicate removal produced 3% more heterozygous SNPs and 12% more heterozygous INDELs before filtering.
  • After filtering, total SNP and INDEL counts were about the same with or without PCR duplicates.The larger pre-filter increase for INDELs suggests greater susceptibility to PCR artifacts than SNPs.

3.6 Effect of the reference genome

Reference-genome completeness strongly affects heterozygous calls: decoy sequences reduce mismapping, while GRCh38 resolves some calls but introduces others of unclear origin.

  • The hs37d5 reference adds 35.4Mb of decoy sequences likely missing from the GRCh37 primary assembly.These sequences are intended to attract reads that would otherwise be mismapped.
  • Homozygous non-low-complexity SNP counts were similar across GRCh37, hs37d5, and GRCh38, but heterozygous call counts differed substantially.The reported homozygous counts were 2.408, 2.405, and 2.412 million, respectively.
  • GRCh37 produced twice as many heterozygotes as hs37d5, indicating that decoy sequences reduced mismapped-read artifacts in chromosomal regions.The authors conclude hs37d5 and GRCh38 are more complete than GRCh37.
  • GRCh38 resolved 39.8k hs37d5 heterozygotes but retained 36.8k GRCh37-only heterozygotes and added 24.6k additional autosomal heterozygotes of unclear origin.These changes show that reference updates can both remove and introduce apparent heterozygous calls.

4 DISCUSSIONS AND CONCLUSIONS

Using haploid CHM1 data, the study identifies major variant-calling artifacts and evaluates practical filtering and consensus strategies. Low-complexity regions and systematic errors remain important challenges, while post-filtering can improve accuracy with little sensitivity loss.

  • Haploid CHM1 enables unbiased identification of erroneous heterozygous calls and testing of filters against diploid NA12878.Because heterozygous calls in CHM1 are generally errors, NA12878 serves as a positive control for sensitivity.
  • 10% of raw CHM1 variant calls were heterozygotes, corresponding to roughly one error per 10–15kb.The authors conclude that most CHM1 heterozygous calls are calling errors.
  • Low-complexity filtering was most effective against false heterozygotes, particularly short INDELs.The result highlights erroneous INDEL realignment in low-complexity regions as a continuing challenge.
  • Filtering low-complexity regions may be practical because over 50% of single-sample INDEL calls occur there, whereas only 1.25% of autosomal ClinVar INDELs overlap them.For some applications, the authors suggest dropping or downweighting calls in these difficult regions.
  • Post-filtered callers usually make one wrong call per 100–200kb outside low-complexity regions without significant sensitivity loss.Different call sets generally agree well on post-filtered sites, although many remaining errors are likely systematic.
  • Intersecting raw calls from two distinct pipelines and applying universal filters can remove most errors with little sensitivity loss.The proposed consensus approach exploits strong agreement on post-filtered sites and poor agreement on false positives.
Loading 1404.0929v2…