IMPORTANT: This is the legacy GATK documentation. This information is only valid until Dec 31st 2019. For latest documentation and forum click here

created by shlee

on 2016-06-21

We are streamlining our recommended workflows by removing a preprocessing step.

As announced in the GATK v3.6 highlights, variant calling workflows that use HaplotypeCaller or MuTect2 now omit indel realignment. This change does not apply to workflows that call variants with UnifiedGenotyper or the original MuTect. We still recommend indel realignment for these legacy workflows. Recall that indel realignment uses RealignerTargetCreator and IndelRealigner and comes after duplicate marking and before base quality score recalibration (BQSR).

In light of these changes, let’s take a brisk stroll through the implications for variant detection. In particular, let’s focus on insertion and deletion events (indels).

Why am I blogging about this now?

I recently studied IndelRealigner and HaplotypeCaller, one after the other, to update a tutorial and to give workshop talks in the UK. This juxtaposition has brought into focus how our capabilities are changing with (i) improving sequencing technology and (ii) improving algorithmic methods. For a comprehensive overview of the state of different sequencing technologies, see the 2016 Goodwin et al. Nature Reviews article. The points in this blogpost pertain in particular to short read data from sequencing by synthesis. Longer short reads, deeper coverage and increased library complexity with PCR-free preparations give us higher quality sequence data. The graph in this external blog shows short read lengths varying from 25 bases in 2006 to 250 bases in 2014. Reads longer than 70 basepairs enable the mapping optimizations BWA-MEM offers, i.e. maximal exact matches. Adding to this, sophisticated haplotype-based variant detectors, such as HaplotypeCaller and FreeBayes, improve our ability to resolve variants. In particular, they are enabling better detection of indels.

Obviously, the added benefit from greater read lengths is more than just the proportional increase in bases. This external article recaps how coverage alone cannot compensate for the shortcomings of short read lengths in closing gaps. Longer reads reduce mapping ambiguity for regions with interspersed repetitive sequences, for some paralogous regions and for gapped alignments. In the case of gapped alignment, the aligner may split the read into multiple alignment records for distant or discordant events, resolve the split with an indel for local events and clip segments that it cannot map. All this splitting and clipping carries some burden of aligner artifact that confounds variant detection. For example, in the case of an insertion in the middle of a read, the inserted sequence distracts from global then local mapping.

In this way, short read mapping limits detection of variants that are dependent on gapped alignment. Towards overcoming these limitations, we added new algorithmic solutions to our workflows --local realignment in 2009 and local assembly in 2012. Illustration A describes a locus with a 20 basepair insertion that we cannot call with just aligner mapped reads. Realigning reads using IndelRealigner or assembling reads using HaplotypeCaller allows us to call the insertion. That indel realignment has been a part of pre-processing workflows for seven years and will continue to be a part of workflows still dependent on locus-based callers is a testament to the improvements it brings. And if you feel apprehensive about omitting it from your HaplotypeCaller and MuTect2 workflows, we empathize. These changes are about improving efficiency in the face of incremental returns. If you find substantial changes, then I encourage you to share details with us.

Let’s consider another illustration that showcases HaplotypeCaller's more robust capability. Illustration B shows how IndelRealigner is unable to realign reads around a 30 basepair insertion, despite inclusion of the region as a target and the presence of a read with the insertion. HaplotypeCaller reassembles and calls this insertion with or without indel realignment. In point of fact, for the example data, the two HaplotypeCaller workflows call identical genotypes for all sites. However, Illustration A points to some subtle differences in the variant annotation metrics. I will talk about the origins of these and their implications.

For those interested in an anecdote from MuTect2, which uses HaplotypeCaller’s approach but without ploidy assumptions to accommodate tumor data, see our MuTect2 workshop presentation. Slides 14–15 outline a case where MuTect2 resolves a 120 basepair deletion haplotype for a tumor.

Did indel realignment still have some value-add?

If indel realignment is redundant to HaplotypeCaller’s functions, why were we recommending both until now? Well, although HaplotypeCaller throws out realignment’s local placement of reads, indel realignment still influences outcomes. I’m told these influences should be subtle for high quality data but make a difference for lower quality data. Let’s talk about the influences conceptually and examine some initial differences for our example data snippet. Breaking down the effects into two categories gives us: (1) the effects of increasing MAPQ by ten for realigned reads and (2) the effects of the physical realignment.

(1) The effects of increasing MAPQ by 10 for realigned reads

(1a) Increases in MAPQ add previously ignored reads to the pool of primary alignment records BaseRecalibrator (BQSR) and HaplotypeCaller consider. BQSR gains previously MAPQ0 reads now MAPQ10 and HaplotypeCaller gains previously MAPQ10–19 reads now over the default MAPQ20 threshold of the HCMappingQualityFilter. (1b) HaplotypeCaller considers all indel realignment target regions, and except for regions that HaplotypeCaller cannot resolve cycles for in its graph assembly, e.g. a tandem repeat region, HaplotypeCaller acts on these regions. Indel realigned reads that cross above HaplotypeCaller's MAPQ threshold alter the allelic and total depth counts for a variant site and can affect genotyping and genotype probability and quality scores as well as other annotation metrics. If the added read(s) are informative towards haplotypes in a region with few informative reads, although they do not impact genotyping any more than uninformative reads, they may have a disproportionate impact on metrics that use exclusively informative reads, e.g. sample-level AD & DP (allelic depth and total depth), QD (QUAL by depth where QUAL is the phred-scaled score for the assertion made in ALT) and FS & SOR (Fisher strand bias and strand odds ratio). (1c) Variant sites that include these increased MAPQ reads directly increase MQ and MQRankSum annotation values (root mean square of mapping quality and ranksum test of ALT vs. REF read mapping qualities). (1d) Finally, these altered annotation metrics (QD, FS, SOR, MQ & MQRankSum for the snippet) impact cohort level annotation metrics, the variant quality score recalibration (VQSR) model and the resulting VQSLOD score for the site. These changes to the distribution of the variant scores impact variant site filtering.

For our data snippet, indel realignment adds zero reads to either pool. The lowest MAPQ of reads indel realignment then increases is at MAPQ26 and in fact over 94% of increases are for MAPQ60 reads. The 194 (~0.07% of total primary) indel realigned reads increase cumulative MAPQ by 0.012% and the impact on MQ for the 20 basepair insertion site in Illustration A is small (MQ 59.39 vs. 62.81). This site has plenty of informative reads. Interestingly, indel realignment adds four more reads that HaplotypeCaller considers for this locus and adds a total of five additional reads for two loci for our snippet's entire one megabasepair region. Of the four additional reads, only one is realigned. This brings us to the next section's topic.

(2) The effects of the physical realignment of reads

Indel realignment shifts read positions and this impacts the outcomes of both BQSR and HaplotypeCaller. (2a) Shifted reads and their mismatching positions around indels will change the bases BQSR modeling uses and potentially recalibrates. Indel realignment generally reduces false positive mismatches around indels. If these mismatches contribute to a covariate pattern, and indel realignment removes these mismatches, then indel realignment removes the support towards these patterns and in turn enables BQSR to calibrate using cleaner models. These contributing false mismatches avoid being unfairly penalized. However, BQSR looks for fairly obvious failures stemming from sequencer flowcell technology such as those this external blog describes neatly. Any discernible patterns from misaligned reads around indels will be small compared to that of flowcell failures. This is one area that could use more data on the extent of the impact of removing indel realignment. What I am told is that effects will focus more on SNPs than indels in that the new workflow lowers the confidence in mainly what will become reference matches after HaplotypeCaller realignment. For these instances, at the worst, you may end up calling a false heterozygous SNP instead of homozygous reference.

Additionally, changes to read positions for a locus impact (2b) the genomic window HaplotypeCaller considers for graph assembly and realignment, because the entropy distribution changes, and (2c) the reads the window includes. For Ilustration A's locus, indel realignment tightens the active region size from 208 to 152 basepairs and also adds four additional reads for consideration. Such occurences alter site annotation values and can influence downstream variant calling.

For our data snippet (again, high quality 30X coverage PCR-free 2x150 WGS), the effects of removing indel realignment are subtle. However, for data that deviate towards lower quality and for data that benefit greatly from BQSR, the effects could be more substantial. Given this, we ask researchers who are considering our new workflows to take into account properties of their data and consider the impact omitting indel realignment may have for their analyses. A first step towards this could be to examine the MAPQ distribution of the primary alignments.

Let's wrap up with some wider picture considerations.

I've covered some ground on the influences of removing indel realignment from a workflow. One side note is that called sites from different variant callers can have different allele variants. For example, IndelRealigner and HaplotypeCaller can differ in how they resolve the variation around a given locus. This raises the question to what extent differential representation, based on different underlying algorithmic solutions, contribute to the discordant calls observed between different workflows. If the extent is great, then, given the potential for differential phenotypic classification, the clinical impact can be significant. A recent article describes how the burgeoning movement to use sequencing in patient diagnostics and disease research is hampered by discordant classifications. By the same reasoning, then the ability to physically phase proximal variants, that determines their cis or trans relationship, has direct impact on phenotypic classification. For example, resolving phasing using tools such as ReadBackedPhasing can tell us whether a variant that shifts the coding frame is expressed on the same molecule as another variant that restores the coding frame.

You should care about resolving indels as best you can even if your analyses focus only on SNPs. Correctly resolving indels reduces the false positive mismatches from alignment that accrue in the vicinity of indels. Also, we care about resolving larger indels for the biological insights they will bring. One example of the biological significance of indels lies in the transposable elements that litter our genomes that in transposing create insertion and deletion events. These elements are currently considered structural given their large size. They are a part of our individual variation and have been proposed to drive speciation and evolution and are implicated in disease including cancer (see endogenous retrovirus). The RepeatMasker web resource catalogues our repetitive elements, to which transposable elements belong, that together comprise over half our human genome.

Finally, technological capabilities are changing the definition of an indel. Back yonder, detection of indels was limited to small indels of a couple to fewer than ten basepairs. GATK workflows assigned resolving anything larger to the realm of structural variant detection. Today, lengthening sequencing reads together with improving global and local alignment algorithms are revealing indels reaching tens to low hundreds of basepairs in our workflows.

To see for yourself the points made in Illustrations A & B, download Tutorial#7156’s example data tutorial7156.tar.gz and follow instructions in Blog#7847README.txt.