#+DATE: [2013-05-06 Mon 08:29]
#+BLOG: bcbio
#+POSTID: 419
#+TITLE: Framework for evaluating variant detection methods: comparison of aligners and callers
#+CATEGORY: variation
#+TAGS: bioinformatics, variant, ngs, clinical, alignment
#+STYLE: <style id="org2blog" type="text/css">td {padding: 4px;}</style>
#+OPTIONS: toc:nil

* Variant detection and grading overview

Developing pipelines for detecting variants from high throughput
sequencing data is challenging due to rapidly changing algorithms and
[[gholson-low][relatively low concordance between methods]]. This post will discuss
automated methods providing evaluation of variant calls, enabling
detailed diagnosis of discordant differences between multiple calling
approaches. This allows us to:

- Characterize strengths and weaknesses of alignment, post-alignment
  preparation and calling methods.

- Automatically verify pipeline updates and installations to ensure
  variant calls recover expected variations. This extends the
  [[xprize-val][XPrize validation protocol]] to provide full summary metrics on
  concordance and discordance of variants.

- Make recommendations on best-practice approaches to use in
  sequencing studies requiring either exome or whole genome variant
  calling.

- Identify characteristics of genomic regions more likely to have
  discordant variants which require additional care when making
  biological conclusions based on calls, or lack of calls, in these
  regions.

This evaluation work is part of a larger community effort to better
characterize variant calling methods.
A key component of these evaluations is a well characterized set of
reference variations for the [[na12878][NA12878 human HapMap genome]], provided by
[[giab][NIST's Genome in a Bottle consortium]]. The diagnostic component of
this work supplements emerging tools like
[[gcat][GCAT (Genome Comparison and Analytic Testing)]],
which provides a community platform for comparing
and discussing calling approaches.

I'll show a 12 way comparison between 2 different aligners
([[novoalign][novoalign]] and [[bwa-mem][bwa mem]]), 2 different post-alignment preparation methods
([[gatk-bp][GATK best practices]] and the [[gkno-me][Marth lab's gkno pipeline]]), and 3
different variant callers ([[gatk-ug][GATK UnifiedGenotyper]],
[[gatk-hc][GATK HaplotypeCaller]], and [[freebayes][FreeBayes]]). This allows comparison of
openly available methods (bwa mem, gkno preparation, and FreeBayes)
with those that require licensing (novoalign, GATK's variant
callers). I'll also describe [[bcbio-nextgen][bcbio-nextgen]], the fully automated open-source
pipeline used for variant calling and evaluation, which
allows others to easily bring this methodology into their
own work and extend this analysis.

#+LINK: gholson-low http://genomemedicine.com/content/5/3/28/abstract
#+LINK: gcat http://www.bioplanet.com/gcat/
#+LINK: ensemble http://bcbio.wordpress.com/2013/02/06/an-automated-ensemble-method-for-combining-and-evaluating-genomic-variants-from-multiple-callers/
#+LINK: giab http://www.genomeinabottle.org/
#+LINK: na12878 http://ccr.coriell.org/Sections/Search/Sample_Detail.aspx?Ref=GM12878
#+LINK: xprize-val http://bcbio.wordpress.com/2012/09/17/genomics-x-prize-public-phase-update-variant-classification-and-de-novo-calling/
#+LINK: novoalign http://www.novocraft.com/main/index.php
#+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: 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: gkno-me http://gkno.me/
#+LINK: freebayes https://github.com/ekg/freebayes

* Aligner, post-alignment preparation and variant calling comparison

To compare methods, we called variants on a NA12878 exome dataset
from [[edge][EdgeBio's clinical pipeline]] and assessed them against the NIST Genome in a
Bottle reference material. Discordant positions where the reference
and evaluation variants differ fall into three different categories:

- Extra variants, called in the evaluation data but not in the
  reference. These are potential false positives.

- Missing variants, found in the NA12878 reference but not in the
  evaluation data set. These are potential false negatives. The use
  of high quality reference materials from NIST enables
  identification of genomic regions we fail to call in.

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

To further identify causes of discordance, we subdivide the missing
and extra variants using annotations from the
[[gemini][GEMINI variation framework]]:

- Low coverage: positions with limited read coverage (4 to 9 reads).
- Repetitive: regions identified by [[repeatmasker][RepeatMasker]].
- Error prone: variants falling in
  [[cse][motifs found to induce sequencing errors]].

We subdivide and restrict our comparisons to help identify sources of
differences between methods indistinguishable when looking at total
discordant counts. A critical subdivison is comparing SNPs and indels
separately. With lower total counts of indels but higher error rates,
each variant type needs independent visualization. Secondly, it's
crucial to distinguish between discordance caused by a lack of
coverage, and discordance caused by an actual difference in
variant assessment. We evaluate only in callable regions with
4 or more reads. This low minimum cutoff provides a valuable
evaluation of low coverage regions, which differ the most between
alignment and calling methods.

I'll use this data to provide recommendations for alignment,
post-alignment preparation and variant calling. In addition to these
high level summaries, the full dataset and summary plots available
below providing a starting place for digging further into the data.

#+LINK: edge http://www.edgebio.com/
#+LINK: gemini https://github.com/arq5x/gemini
#+LINK: repeatmasker http://repeatmasker.org/
#+LINK: cse http://www.biomedcentral.com/1471-2105/14/S5/S1

** Aligners

We compared two recently released aligners designed to work with
longer reads coming from new sequencing technologies:
[[novoalign][novoalign (3.00.02)]] and [[bwa-mem][bwa mem (0.7.3a]]). bwa mem identified 1389 additional
concordant SNPs and 145 indels not seen with novoalign. 1024 of
[[bwa-novo-diffs][these missing variants]] are in regions where novoalign does not provide
sufficient coverage for calling. Of those, 92% (941) have low coverage
with less than 10 reads in the bwa alignments. Algorithmic changes impact low
coverage regions more due to the decreased evidence and
susceptibility to crossing calling coverage thresholds, so we need
extra care and consideration of calls in these regions.

Our standard workflow uses novoalign based on its stringency in
[[align-paper][resolving large insertions and deletions]]. These results suggest
equally good results using bwa mem, along with improved processing
times. One caveat to these results is that some of the available
Illumina call data that feeds into NIST's reference genomes comes
from a bwa alignment, so some differences may reflect a bias
towards bwa alignment heuristics. Using non-simulated
reference data sets has the advantage of capturing real biological
and process errors, but requires iterative improvement of the
reference materials to avoid this type of potential algorithmic bias.

#+ATTR_HTML: alt="Comparison of concordant variants by aligner type" width="700"
[[file:grading-summary-prep-alignerdiff.png]]

#+grading-summary-prep-alignerdiff.png https://raw.github.com/chapmanb/bcbb/master/posts/calling_pipeline_compare/grading-summary-prep-alignerdiff.png
#+LINK: align-paper http://f1000research.com/articles/1-2/v2#f3

** Post-alignment preparation and quality score recalibration

We compared two methods of quality recalibration:

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

- [[gkno-me][The Marth Lab's gkno realignment pipeline]]: This performs de-duplication
  with [[samtools][samtools rmdup]] and realignment around indels using [[ogap][ogap]]. All
  commands in this pipeline work on streaming input, avoiding disk IO
  penalties by using unix pipes. This piped approach improves
  scaling on large numbers of whole genome samples. Notably, our
  implementation of the pipeline does not use a base quality score
  recalibration step.

GATK best practice pipelines offer an advantage over the gkno-only
pipeline primarily because of improvements in SNP calling from base
quality recalibration. Specifically we lose ~1% (824 / 77158) of
called variations. These fall into the discordant missing "other"
category, so we cannot explain them by metrics such as coverage or
genome difficulty. The simplest explanation is that initial poor
quality calculations in those regions result in callers missing those
variants. Base quality recalibration helps recover them. These results match
[[brendan-qual][Brendan O'Fallon's recent analysis of base quality score recalibration]].

This places a practical number on the lost variants when avoiding
recalibration either due to scaling or GATK licensing concerns. Some
other options for recalibration include
[[novoalign-qual][Novoalign's Quality Recalibration]] and
[[bamtuil][University of Michigan's BamUtil recab]], although we've not yet tested
either in depth as potential supplements to improve calling in
non-GATK pipelines.

#+ATTR_HTML: alt="Comparison of concordant variants by post-alignment prep method" width="700"
[[file:grading-summary-prep-bamprepdiff.png]]

#+grading-summary-prep-bamprepdiff.png https://raw.github.com/chapmanb/bcbb/master/posts/calling_pipeline_compare/grading-summary-prep-bamprepdiff.png
#+LINK: brendan-qual http://basecallbio.wordpress.com/2013/04/23/base-quality-score-rebinning/
#+LINK: samtools http://samtools.sourceforge.net/
#+LINK: ogap https://github.com/ekg/ogap
#+LINK: picard-md http://picard.sourceforge.net/command-line-overview.shtml#MarkDuplicates
#+LINK: novalign-qual http://novocraft.com/wiki/tiki-index.php?page=Quality+Calibration
#+LINK: bamutil http://genome.sph.umich.edu/wiki/BamUtil

** Variant callers

For this comparison, we used three general purpose callers that handle
SNPs and small indels, all of which have updated versions since
[[ensemble][our last comparison]]:

- [[freebayes][FreeBayes (0.9.9 296a0fa)]]: A haplotype-based Bayesian caller from
  the Marth Lab, with filtering on quality score and read depth.

- [[gatk-ug][GATK UnifiedGenotyper (2.4-9)]]: GATK's widely used Bayesian caller,
  using filtering recommendations for exome experiments from
  [[gatk-bp][GATK's best practices]].

- [[gatk-hc][GATK HaplotypeCaller (2.4-9)]]: GATK's more recently developed
  haplotype caller which provides local assembly around variant
  regions, using filtering recommendations for exomes from
  [[gatk-bp][GATK's best practices]].

Adjusting variant calling methods has the biggest impact on the final
set of calls. Called SNPs differ by 4577 between the three compared
approaches, in comparison with aligner and post-alignment preparation
changes which resulted in a maximum difference of 1389 calls. This
suggests that experimenting with variant calling approaches currently
provides the most leverage to improve calls.

A majority of the SNP concordance differences between the three calling
methods are in low coverage regions with between 4 and 9
reads. GATK UnifiedGenotyper performs the best in detecting SNPs in
these low coverage regions. FreeBayes and GATK HaplotypeCaller both
call more conservatively in these regions, generating more potential
false negatives. FreeBayes had the fewest heterozygote/homozygote
discrimination differences of the three callers.

For indels, FreeBayes and HaplotypeCaller both provide improved
sensitivity compared to UnifiedGenotyper, with HaplotypeCaller
identifying the most, especially in low coverage regions. In contrast
to the SNP calling results, FreeBayes has more calls that match the
expected indel but differ in whether a call is a heterozygote or
homozygote.

#+ATTR_HTML: alt="Comparison of concordant variants by calling method" width="800"
[[file:grading-summary-prep-callerdiff.png]]

No one caller outperformed the others on all subsets of the data. GATK
UnifiedGenotyper performs best on SNPs but is less sensitive in
resolving indels. GATK HaplotypeCaller identifies the most indels, but
is more conservative than the other callers on SNPs. FreeBayes
provides intermediate sensitivity and specificity between the two for
both SNPs and indels. A combined UnifiedGenotyper and HaplotypeCaller
pipeline for SNPs and indels, respectively, would provide the best
overall calling metrics based on this set of comparisons.

Low coverage regions are the key area of difference between
callers. Coupled with the alignment results and investigation of
variant changes resulting from [[bcbio-qb][quality score binning]], this suggests we should
be more critical in assessing both calls and coverage in these
regions. Assessing coverage and potential false negatives is
especially critical since we lack good tools to summarize and
prioritize genomic regions that are potentially missed during
sequencing. This also emphasizes the role of population-based calling
to help resolve low coverage regions, since callers can use evidence
from multiple samples to better estimate the likelihoods of low
coverage calls.

#+LINK: bcbio-qb http://bcbio.wordpress.com/2013/02/13/the-influence-of-reduced-resolution-quality-scores-on-alignment-and-variant-calling/
#+grading-summary-prep-callerdiff.png https://raw.github.com/chapmanb/bcbb/master/posts/calling_pipeline_compare/grading-summary-prep-callerdiff.png

* Automated calling and grading pipeline

Method comparisons become dated quickly due to the continuous
improvement in aligners and variant callers. While these
recommendations are useful now, in 6 months there will be new releases
with improved approaches. This rapid development cycle creates
challenges for biologists hoping to derive meaning from variant
results: do you stay locked on software versions whose trade offs you
understand, or do you attempt to stay current and handle re-verifying
results with every new release?

Our goal is to provide a community developed pipeline and comparison
framework that ameliorates this continuous struggle to re-verify. The
analysis done here is fully automated as part of the [[bcbio-nextgen][bcbio-nextgen]]
analysis framework. This framework code provides full exposure
and revision tracking of all parameters used in analyses. For example,
the ngsalign module contains the command lines used for
[[cl-bwa][bwa mem]] and [[cl-novoalign][novoalign]], as well as all other tools.

To install the pipeline, third-party software and required data files:

#+BEGIN_EXAMPLE
    wget https://raw.github.com/chapmanb/bcbio-nextgen/master/scripts/bcbio_nextgen_install.py
    python bcbio_nextgen_install.py /usr/local /usr/local/share/bcbio-nextgen
#+END_EXAMPLE
#+begin_html
 <br />
#+end_html

The installer bootstraps all installation on a bare machine using
[[cbl][the CloudBioLinux framework]]. More details and options are available
in the [[install-docs][installation documentation]].

To re-run this analysis, retrieve the input data files and configuration
as described in the [[example-docs][bcbio-nextgen example documentation]] with:

#+BEGIN_EXAMPLE
    $ mkdir config && cd config
    $ wget https://raw.github.com/chapmanb/bcbio-nextgen/master/config/\
       examples/NA12878-exome-methodcmp.yaml
    $ cd .. && mkdir input && cd input
    $ wget https://dm.genomespace.org/datamanager/file/Home/EdgeBio/\
       CLIA_Examples/NA12878-NGv3-LAB1360-A/NA12878-NGv3-LAB1360-A_1.fastq.gz
    $ wget https://dm.genomespace.org/datamanager/file/Home/EdgeBio/\
       CLIA_Examples/NA12878-NGv3-LAB1360-A/NA12878-NGv3-LAB1360-A_2.fastq.gz
    $ wget https://s3.amazonaws.com/bcbio_nextgen/NA12878-nist-v2_13-NGv3-pass.vcf.gz
    $ wget https://s3.amazonaws.com/bcbio_nextgen/NA12878-nist-v2_13-NGv3-regions.bed.gz
    $ gunzip NA12878-nist-*.gz
    $ wget https://s3.amazonaws.com/bcbio_nextgen/NGv3.bed.gz
    $ gunzip NGv3.bed.gz
#+END_EXAMPLE
#+begin_html
 <br />
#+end_html

Then run the analysis, distributed on 8 local cores, with:

#+BEGIN_EXAMPLE
    $ mkdir work && cd work
    $ bcbio_nextgen.py bcbio_system.yaml ../input ../config/NA12878-exome-methodcmp.yaml -n 8
#+END_EXAMPLE
#+begin_html
 <br />
#+end_html

The [[parallel-docs][bcbio-nextgen documentation]] describes how to
parallelize processing over multiple machines using cluster
schedulers (LSF, SGE, Torque).

The pipeline and comparison framework are open-source and
configurable for multiple aligners, preparation methods and callers.
We invite anyone interested in this work to provide feedback and
contributions.

#+LINK: bcbio-nextgen https://github.com/chapmanb/bcbio-nextgen
#+LINK: cbl http://cloudbiolinux.org
#+LINK: install-docs https://bcbio-nextgen.readthedocs.org/en/latest/contents/installation.html
#+LINK: parallel-docs https://bcbio-nextgen.readthedocs.org/en/latest/contents/parallel.html
#+LINK: example-docs https://bcbio-nextgen.readthedocs.org/en/latest/contents/testing.html#example-pipelines
#+LINK: cl-novoalign https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/ngsalign/novoalign.py#L44
#+LINK: cl-bwa https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/ngsalign/bwa.py#L38

* Full data sets

We extracted the conclusions for alignment, post-alignment
preparation and variant calling from analysis of the full
dataset. The visualizations for the full data are not as pretty but we
make them available for anyone interested in digging deeper:

- [[summary-csv][Summary CSV of comparisons]] split by methods and
  concordance/discordance types, easily importable into [[r][R]] or [[pandas][pandas]]
  for further analysis.
- [[plot-code][Code for preparing and plotting results]]
- Full comparisons of all 12 methods, stratified by concordance and
  discordance: [[full-snp][SNPs]] and [[full-indel][indels]]
- Boxplots of differences between alignment methods: [[alignerdiff-snp][SNPs]] and [[alignerdiff-indel][indels]]
- Boxplots of differences between post-alignment preparation methods:
  [[bamprepdiff-snp][SNPs]] and [[bamprepdiff-indel][indels]]
- Boxplots of differences between variant calling methods: [[callerdiff-snp][SNPs]] and [[callerdiff-indel][indels]]

The comparison variant calls are also useful for pinpointing algorithmic
differences between methods. Some useful subsets of variants:

- Concordant variants called by bwa and not novoalign, where novoalign did
  not have sufficient coverage in the region. These are calls where
  either novoalign fails to map some reads, or bwa maps too aggressively:
  [[bwa-novo-diffs][VCF of bwa calls with low or no coverage in novoalign.]]

- Discordant variants called consistently by multiple calling methods.
  These are potential errors in the reference material, or
  consistently problematic calling regions for multiple
  algorithms. Of the 9004 shared discordants, the majority are
  potential false negatives not seen in the evaluation calls (7152; 79%). Another large portion is
  heterozygote/homozygote differences, which make up 1627 calls (18%).
  6652 (74%) of the differences have low coverage in the exome
  evaluation, again reflecting the difficulties in calling in these regions.
  The [[shared-discordants][VCF of discordants found in 2 or more callers]] contains these
  calls, with a 'GradeCat' INFO tag specifying the discordance category.

We encourage reanalysis and welcome suggestions for improving the presentation
and conclusions in this post.

#+LINK: r http://cran.r-project.org/
#+LINK: pandas http://pandas.pydata.org/
#+LINK: summary-csv https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep.csv
#+LINK: plot-code https://github.com/chapmanb/bcbb/tree/master/validation
#+LINK: full-snp https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-SNP.pdf
#+LINK: full-indel https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-SNP.pdf
#+LINK: alignerdiff-snp https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-aligner-SNP.pdf
#+LINK: alignerdiff-indel https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-aligner-Indel.pdf
#+LINK: bamprepdiff-snp https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-bamprep-SNP.pdf
#+LINK: bamprepdiff-indel https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-bamprep-Indel.pdf
#+LINK: callerdiff-snp https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-caller-SNP.pdf
#+LINK: callerdiff-indel https://github.com/chapmanb/bcbb/raw/master/posts/calling_pipeline_compare/grading-summary-prep-caller-Indel.pdf
#+LINK: bwa-novo-diffs https://dm.genomespace.org/datamanager/file/Home/chapmanb/validation/NA12878_cmp/NA12878-cmp-in-bwa-not-novoalign.vcf
#+LINK: shared-discordants https://dm.genomespace.org/datamanager/file/Home/chapmanb/validation/NA12878_cmp/NA12878-cmp-shared-discordants.vcf
