HaplotypeCaller in a nutshell

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 GATK_Team

on 2017-12-28

This document outlines the basic operation of the HaplotypeCaller run in its default mode on a single sample, and does not cover the additional processing and calculations done when it is run in "GVCF mode" (with -ERC GVCF or -ERC BP_RESOLUTION) or when it is run on multiple samples. For more details and discussion of the GVCF workflow, see the Best Practices documentation on germline short variant discovery as well as the HaplotypeCaller manuscript on bioRxiv.

Overview

The core operations performed by HaplotypeCaller can be grouped into these major steps:

1. Define active regions. The program determines which regions of the genome it needs to operate on, based on the presence of significant evidence for variation.

2. Determine haplotypes by re-assembly of the active region. For each ActiveRegion, the program builds a De Bruijn-like graph to reassemble the ActiveRegion and identifies what are the possible haplotypes present in the data. The program then realigns each haplotype against the reference haplotype using the Smith-Waterman algorithm in order to identify potentially variant sites.

3. Determine likelihoods of the haplotypes given the read data. For each ActiveRegion, the program performs a pairwise alignment of each read against each haplotype using the PairHMM algorithm. This produces a matrix of likelihoods of haplotypes given the read data. These likelihoods are then marginalized to obtain the likelihoods of alleles per read for each potentially variant site.

4. Assign sample genotypes. For each potentially variant site, the program applies Bayes’ rule, using the likelihoods of alleles given the read data to calculate the posterior likelihoods of each genotype per sample given the read data observed for that sample. The most likely genotype is then assigned to the sample.

1. Define active regions

In this first step, the program traverses the sequencing data to identify regions of the genomes in which the samples being analyzed show substantial evidence of variation relative to the reference. The resulting areas are defined as “active regions”, and will be passed on to the next step. Areas that do not show any variation beyond the expected levels of background noise will be skipped in the next step. This aims to accelerate the analysis by not wasting time performing reassembly on regions that are identical to the reference anyway.

To define these active regions, the program operates in three phases. First, it computes an activity score for each individual genome position, yielding the raw activity profile, which is a wave function of activity per position. Then, it applies a smoothing algorithm to the raw profile, which is essentially a sort of averaging process, to yield the actual activity profile. Finally, it identifies local maxima where the activity profile curve rises above the preset activity threshold, and defines appropriate intervals to encompass the active profile within the preset size constraints. For more details on how the activity profile is computed and processed, as well as what options are available to modify the active region parameters, please see this article.

Once this process is complete, the program applies a few post-processing steps to finalize the the active regions (see detailed doc above). The final output of this process is a list of intervals corresponding to the active regions which will be processed in the next step.

2. Determine haplotypes by local assembly of the active region.

The goal of this step is to reconstruct the possible sequences of the real physical segments of DNA present in the original sample organism. To do this, the program goes through each active region and uses the input reads that mapped to that region to construct complete sequences covering its entire length, which are called haplotypes. This process will typically generate several different possible haplotypes for each active region due to:

• real diversity on polyploid (including CNV) or multi-sample data
• possible allele combinations between variant sites that are not totally linked within the active region
• sequencing and mapping errors

In order to generate a list of possible haplotypes, the program first builds an assembly graph for the active region using the reference sequence as a template. Then, it takes each read in turn and attempts to match it to a segment of the graph. Whenever portions of a read do not match the local graph, the program adds new nodes to the graph to account for the mismatches. After this process has been repeated with many reads, it typically yields a complex graph with many possible paths. However, because the program keeps track of how many reads support each path segment, we can select only the most likely (well-supported) paths. These likely paths are then used to build the haplotype sequences which will be used for scoring and genotyping in the next step.

The assembly and haplotype determination procedure is described in full detail in this method article.

Once the haplotypes have been determined, each one is realigned against the original reference sequence in order to identify potentially variant sites. This produces the set of sites that will be processed in the next step. A subset of these sites will eventually be emitted as variant calls to the output VCF.

3. Evaluating the evidence for haplotypes and variant alleles

Now that we have all these candidate haplotypes, we need to evaluate how much evidence there is in the data to support each one of them. So the program takes each individual read and aligns it against each haplotype in turn (including the reference haplotype) using the PairHMM algorithm, which takes into account the information we have about the quality of the data (i.e. the base quality scores and indel quality scores). This outputs a score for each read-haplotype pairing, expressing the likelihood of observing that read given that haplotype.

Those scores are then used to calculate out how much evidence there is for individual alleles at the candidate sites that were identified in the previous step. The process is called marginalization over alleles and produces the actual numbers that will finally be used to assign a genotype to the sample in the next step.

For further details on the pairHMM output and the marginalization process, see this document.

4. Assigning per-sample genotypes

The previous step produced a table of per-read allele likelihoods for each candidate variant site under consideration. Now, all that remains to do is to evaluate those likelihoods in aggregate to determine what is the most likely genotype of the sample at each site. This is done by applying Bayes' theorem to calculate the likelihoods of each possible genotype, and selecting the most likely. This produces a genotype call as well as the calculation of various metrics that will be annotated in the output VCF if a variant call is emitted.

For further details on the genotyping calculations, see this document.

This concludes the overview of how HaplotypeCaller works.

Updated on 2018-02-23

From blueskypy on 2018-09-13

Seems all the documents of HC does not distinguish the difference between “probability” and “likelihood”, and uses “likelihood” often in the place of “probability”.

P(Gi|D) is the (posterior) probability of the genotype i given the observed reads

P(D|Gi) is the likelihood of Gi given the observed reads.

Not distinguishing the two terms makes reading the step 4 very confusing!

From shlee on 2018-09-24

Hey @blueskypy,

I’m just helping out our new forum support person today and I’m noticing that you’ve posted four similar questions across different threads for the last few days. Would it be possible for you to consolidate all of your questions into a [single new post](https://gatkforums.broadinstitute.org/gatk/post/question/gatk-support-forum)?

Here are your four threads:

- https://gatkforums.broadinstitute.org/gatk/discussion/comment/52389#Comment_52389

- https://gatkforums.broadinstitute.org/gatk/discussion/comment/52386#Comment_52386

- https://gatkforums.broadinstitute.org/gatk/discussion/comment/52368#Comment_52368

- https://gatkforums.broadinstitute.org/gatk/discussion/comment/52183#Comment_52183

Going forward, it would help us tremendously if you could not ask duplicate questions and to group your related questions into the same thread. Thank you!

From Junyan on 2018-11-12

Hi, I have a simple question about haplotypercaller.

When run HaplotypeCaller with parameter -I input.bam and -L interval_list, these core operations (like define active region, local assembly, determine haplotype …) performed by HaplotypeCaller are based on genome region in bam or Interval list provided? The Interval list is used in the whole Haplotype process or only used to filter variants?

Thanks a lot.