#+DATE: [2013-10-21 Mon 06:35]
#+BLOG: bcbio
#+POSTID: 540
#+TITLE: Updated comparison of variant detection methods: Ensemble, FreeBayes and minimal BAM preparation pipelines
#+CATEGORY: variation
#+TAGS: bioinformatics, variant, ngs, clinical
#+OPTIONS: toc:nil num:nil

* Variant evaluation overview

I previously discussed our approach for [[eval-variant][evaluating variant detection methods]]
using a [[giab-paper][highly confident set of reference calls]] provided by
[[giab][NIST's Genome in a Bottle consortium]] for the [[na12878][NA12878 human HapMap genome]],
In this post, I'll update those conclusions based on recent improvements
in GATK and FreeBayes.

The comparisons use [[bcbio-nextgen][bcbio-nextgen]], an automated open-source
pipeline for variant calling and evaluation that identifies concordant
and discordant variants with the [[xprize-val][XPrize validation protocol]].
By having an automated validation workflow attached to a
regularly updated, community developed, variant calling pipeline, we
can actively track progress of variant callers and provide updates as
algorithms improve.

Since the initial post, There have been two new GATK releases of
[[gatk-ug][UnifiedGenotyper]] and [[gatk-hc][HaplotypeCaller]], as well as multiple improvements
to [[freebayes][FreeBayes]]. Additionally we've enchanced our [[ensemble][ensemble calling method]],
which combines inputs from multiple callers into a single
final set of calls, to better handle comparisons with inputs from
three callers.

The goal of this post is to re-evaluate these variant detection
approaches and provide an updated set of recommendations:

- FreeBayes detects more concordant SNPs and indels compared to GATK
  approaches, including GATK's HaplotypeCaller method.

- Post-alignment BAM processing steps like base quality recalibration and
  realignment have little impact on the quality of variant calls with
  variant callers that perform local realignment, including FreeBayes
  and GATK HaplotypeCaller.

- The Ensemble calling method provides the best variant detection by
  combining inputs from GATK UnifiedGenotyper, HaplotypeCaller and
  FreeBayes.

Avoiding the post-alignment BAM recalibration and realignment steps
allows us to save significant time and pipeline complexity.
Combined with the improvements in FreeBayes, this enables
a variant calling pipeline that can be freely used for academic,
clinical and commercial work with equal quality variant calls compared
to current GATK best-practice approaches.

#+LINK: eval-variant http://bcbio.wordpress.com/2013/05/06/framework-for-evaluating-variant-detection-methods-comparison-of-aligners-and-callers/
#+LINK: na12878 http://ccr.coriell.org/Sections/Search/Sample_Detail.aspx?Ref=GM12878
#+LINK: giab-paper http://arxiv.org/abs/1307.4661
#+LINK: giab http://www.genomeinabottle.org/
#+LINK: xprize-val http://bcbio.wordpress.com/2012/09/17/genomics-x-prize-public-phase-update-variant-classification-and-de-novo-calling/
#+LINK: freebayes https://github.com/ekg/freebayes
#+LINK: gatk-ug http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_genotyper_UnifiedGenotyper.html
#+LINK: gatk-hc http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_haplotypecaller_HaplotypeCaller.html
#+LINK: ensemble http://bcbio.wordpress.com/2013/02/06/an-automated-ensemble-method-for-combining-and-evaluating-genomic-variants-from-multiple-callers/
#+LINK: bcbio-nextgen https://github.com/chapmanb/bcbio-nextgen

* Calling and evaluation methods

We called variants on a NA12878 exome dataset
from [[edge][EdgeBio's clinical pipeline]] and assessed them against the NIST's Genome in a
Bottle reference material. [[comparison-do][Full instructions for replicating the analysis]]
and [[bcbio-install][installing the pipeline]] are available from the bcbio-nextgen documentation site.
Following alignment with [[bwa-mem][bwa-mem (0.7.5a)]], we post-processed the BAM
files with two methods:

- [[gatk-bp][GATK's best practices (2.7-2)]]: This involves de-duplication with
  [[picard-md][Picard MarkDuplicates]], [[gatk-bqsr][GATK base quality score recalibration]] and
  [[gatk-realign][GATK realignment around indels]].

- Minimal post-processing, with de-duplication using
  [[samtools][samtools rmdup]] and no realignment or recalibration.

We then called variants with three general purpose callers:

- [[freebayes][FreeBayes (v0.9.9.2-18)]]: A haplotype-based Bayesian caller from
  the Marth Lab. We filter calls with a hard filter based on depth,
  quality and strand bias.

- [[gatk-ug][GATK UnifiedGenotyper (2.7-2)]]: GATK's widely used Bayesian caller.
  Since this is single sample exome data, we filter calls using
  [[gatk-hardfilter][GATK's recommended hard filters]], instead of
  [[broad-vqsr][Variant Quality Score Recalibration (VQSR)]].

- [[gatk-hc][GATK HaplotypeCaller (2.7-2)]]: GATK's more recently developed
  haplotype caller which provides local assembly around variant
  regions. We also filtered these calls using recommended hard filters.

Finally, we evaluated the calls from each combination of variant
caller and BAM post-alignment preparation method using the
[[bcbio.variation][bcbio.variation]] framework. This provides a summary identifying
concordant and discordant variants, separating SNPs and indels since
they have different error profiles. Additionally it classifies
discordant variants. where the reference material and evaluation
variants differ, into three categories:

- Extra variants, called in the evaluation data but not in the
  reference. These are potential false positives or missing calls from
  the reference materials.

- Missing variants, found in the NA12878 reference but not in the
  evaluation data set. These are potential false negatives.

- Shared variants, called in both the evaluation and reference but
  differently represented. This results from allele differences, such as
  heterozygote versus homozygote calls, or variant identification
  differences, such as indel start and end coordinates.

#+LINK: edge http://www.edgebio.com/
#+LINK: bwa-mem http://bio-bwa.sourceforge.net/
#+LINK: gatk-bp http://gatkforums.broadinstitute.org/discussion/1186/best-practice-variant-detection-with-the-gatk-v4-for-release-2-0
#+LINK: comparison-do https://bcbio-nextgen.readthedocs.org/en/latest/contents/testing.html#exome-with-validation-against-reference-materials
#+LINK: samtools http://samtools.sourceforge.net/
#+LINK: picard-md http://picard.sourceforge.net/command-line-overview.shtml#MarkDuplicates
#+LINK: bcbio.variation https://github.com/chapmanb/bcbio.variation
#+LINK: broad-vqsr http://gatkforums.broadinstitute.org/discussion/39/variant-quality-score-recalibration-vqsr
#+LINK: gatk-hardfilter http://gatkforums.broadinstitute.org/discussion/2806/howto-apply-hard-filters-to-a-call-set
#+LINK: bcbio-install https://bcbio-nextgen.readthedocs.org/en/latest/contents/installation.html#automated

* Variant caller comparison

Using this framework, we compared the 3 variant callers and combined
ensemble method:

- FreeBayes outperforms the GATK callers on both SNP and indel
  calling. The most recent versions of FreeBayes have improved
  sensitivity and specificity which puts them on par with GATK
  HaplotypeCaller. One area where FreeBayes performs better is in
  correctly resolving heterozygote/homozygote calls, reflected in the
  lower number of discordant shared variants.

- GATK HaplotypeCaller is all around better than the UnifiedGenotyper.
  In the previous comparison, we found UnifiedGenotyper performed
  better on SNPs and HaplotypeCaller better on indels, but the recent
  improvements in GATK 2.7 have resolved the difference in SNP
  calling. If using a GATK pipeline, UnifiedGenotyper lags behind the
  realigning callers in resolving indels, and I'd recommend using
  HaplotypeCaller. This mirrors the GATK team's current recommendations.

- The ensemble calling approach provides the best overall resolution
  of both SNPs and indels. The one area where it lags slightly behind
  is in identification of homozygote/heterozygote calls, especially in
  indels. This is due to positions where HaplotypeCaller and FreeBayes
  both call variants but differ on whether it is a heterozygote or
  homozygote, reflected as higher discordant shared counts.

#+BEGIN_HTML
<a href="http://i.imgur.com/qz4Maf6.png">
  <img src="http://i.imgur.com/qz4Maf6.png" width="700"
       alt="Comparison of variant callers, GATK best practice preparation">
</a>
#+END_HTML

In addition to calling sensitivity and specificity, an additional
factor to consider is the required processing time. Rough benchmarks
on family-based calling of whole genome sequencing data indicate that
HaplotypeCaller is roughly 7x slower than UnifiedGenotyper and
FreeBayes is 2x slower. On multiple 30x whole genome samples, our
experience is that calling can range from 10 hours for GATK
UnifiedGenotyper to 70 hours for HaplotypeCallers. Ensemble calling
requires running all three callers plus combining into a final call
set, and for family-based whole genome samples can add another
100 hours of processing time. These estimates fluctuate greatly
depending on the compute infrastructure and presence of longer
difficult genomic regions with deeper coverage, but give some
estimates of timing considerations.

* Post-alignment BAM preparation comparison

Given the improved accuracy of local realignment haplotype-based
callers like FreeBayes and HaplotypeCaller, we explored the accuracy
cost of removing the post-alignment BAM processing steps. The
recommended GATK best-practice is to follow up alignment with
identification of duplicate reads, followed by
[[gatk-bqsr][base quality score recalibration]] and [[gatk-realign][realignment around indels.]]
Based on [[bcbio-scale][whole genome benchmarking work]], these steps can take as long
as the initial alignment and scale poorly due to the high IO costs of
manipulating large BAM files. For multiple 30x whole genome samples
running on 16 cores per sample, this can account for 12 to 16 hours of
processing time.

To compare the quality impact of avoiding recalibration and
realignment, we performed the identical alignment and variant calling
steps as above, but did minimal post-alignment BAM preparation.
Following alignment, the only step performed was deduplication using
[[samtools][samtools rmdup]]. Unlike Picard MarkDuplicates, samtools rmdup
handles piped streaming input to avoid IO penalties. This is
at the [[rmdup-v-markdup][cost of not handling some edge cases]]. Longer term, we'd like to
explore [[biobambam][biobambam's markduplicates2]], which implements a more efficient
streaming version of the Picard MarkDuplicates algorithm.

Suprisingly, skipping base recalibration and indel realignment had
almost no impact on the quality of resulting variant calls:

#+BEGIN_HTML
<a href="http://i.imgur.com/w8g0HCv.png]">
  <img src="http://i.imgur.com/w8g0HCv.png]" width="700"
       alt="Comparison of variant callers, minimal post-alignment preparation">
</a>
#+END_HTML

While GATK UnifiedGenotyper suffers during indel calling without
recalibration and realignment, both HaplotypeCaller and FreeBayes
perform as good or better without these steps. This allows us to save
on processing time and complexity without sacrificing call quality
when using a haplotype aware realigning caller.

#+LINK: gatk-bqsr http://gatk.vanillaforums.com/discussion/44/base-quality-score-recalibration-bqsr
#+LINK: gatk-realign http://gatk.vanillaforums.com/discussion/38/local-realignment-around-indels
#+LINK: bcbio-scale http://bcbio.wordpress.com/2013/05/22/scaling-variant-detection-pipelines-for-whole-genome-sequencing-analysis/
#+LINK: biobambam https://github.com/gt1/biobambam
#+LINK: rmdup-v-markdup http://www.biostars.org/p/3917/#3985

* Caveats and conclusions

Taken together, the improvements in FreeBayes and ability to avoid
post-alignment BAM processing allow use of a commercially unrestricted
GATK-free pipeline with equal quality to current GATK best practices.
Adding in GATK's two callers plus our ensemble combining method
provides the most accurate overall calls, at the cost of
additional processing time.

It's also important to consider potential drawbacks of this analysis
as we continue to design future evaluations. The comparison is in
exome regions for single sample variant calling. In future work it
would be helpful to have population or family based inputs. We'd also
like to prepare test datasets that focus specifically on evaluating
the quality of calls in more difficult repetitive regions within the
whole genome. Using populations or whole genomes would also allow use
of GATK's Variant Quality Score Recalibration as part of the pipeline,
which could provide improved filtering compared to the hard-filtering
approach used here.

Another consideration is that the reference callset prepared by the
Genome in a Bottle consortium makes extensive use of GATK tools
during preparation. Evaluation of the reference materials with
FreeBayes and other callers can help reduce potential GATK-specific
biases when continuing to develop reliable reference materials.

All of these pipelines are freely available, open-source, community
developed projects and we welcome feedback and contributors. By
integrating validation into a scalable analysis pipeline, we hope to
build a community interested in widely accessible calling pipelines
coupled with well-evaluated reference datasets and methods.
