#+BLOG: bcbio
#+POSTID: 748
#+DATE: [2015-03-05 Thu 06:13]
#+TITLE: Validating multiple cancer variant callers and prioritization in tumor-only samples
#+CATEGORY: cancer
#+TAGS: bioinformatics, variant, ngs, validation, cancer, somatic, mutect, freebayes, vardict, scalpel
#+OPTIONS: toc:nil num:nil

* Overview

The post discusses work validating multiple cancer variant callers in
[[bcbio][bcbio-nextgen]] using a synthetic reference call set from
[[dream][the ICGC-TCGA DREAM challenge]].  We've previously [[jointval][validated germline variant calling methods]],
but cancer calling is additionally challenging. Tumor samples have mixed cellularity
due to contaminating normal sample, and consist of multiple sub-clones with
different somatic variations. Low-frequency sub-clonal variations can be
critical to understand disease progression but are more difficult to detect with
high sensitivity and precision.

Publicly available whole genome truth sets like the
[[giab][NA12878 Genome in a Bottle reference materials]] don't currently exist for cancer
calling, but other groups have been working to prepare standards for use in evaluating
callers. The [[dream][DREAM challenge]] provides a set of [[dream_synthetic][synthetic datasets]] that include cellularity and multiple
sub-clones. There is also on ongoing DREAM contest with real, non-simulated
data. In addition, the [[icgc][ICGC]] has done extensive work on
[[icgc_calling][assessing the challenges in cancer variant calling]] and produced a detailed
comparison of [[icgc_pipelines][multiple variant calling pipelines]]. Finally, [[bina][Bina]] recently
released a [[bina_varsim][simulation framework called varsim]] that they used to
[[bina_varsim_cancer][compare multiple callers]] as part of their [[bina_cancer][cancer benchmarking]]. I'm excited
about all this community benchmarking and hope this work can contribute to the
goal of having fully consented real patient reference materials, leading to
a good understanding of variant caller tradeoffs.

In this work, we [[dream_about][evaluated cancer tumor/normal variant calling]] with
[[dream_synthetic][synthetic dataset 3]] from the DREAM challenge, using multiple approaches to detect SNPs,
indels and structural variants. A majority voting ensemble approach combines inputs from
multiple callers into a final callset with good sensitivity and precision. We
also provide a prioritization method to enrich for somatic mutations in tumor
samples without matched normals.

Cancer variant calling in [[bcbio][bcbio]] is due to contributions and support from
multiple members of the community:

- [[miika][Miika Ahdesmaki]] and [[justin][Justin Johnson]] at [[az][AstraZeneca]] collaborated with our group
  on integrating and evaluating multiple variant callers. Their financial
  support helped fund our time on this comparison, and their scientific
  contributions improved SNP and indel calling.

- [[luca][Luca Beltrame]] integrated a number of cancer variant callers into bcbio and
  wrote the initial framework for somatic calling.

- [[lorena][Lorena Pantano]] integrated variant callers and performed benchmarking
  work. Members of our group, the [[hbc][Harvard Chan school Bioinformatics core]],
  regularly contribute to bcbio development and validation work.

- The [[wwcrc][Wolfson Wohl Cancer Research Centre]] supported our work on validation of
  somatic callers.

- James Cuff, Paul Edmon and the team at [[https://rc.fas.harvard.edu/][Harvard FAS research computing]]
  provided compute infrastructure and support that enabled the large number
  of benchmarking runs it takes to get good quality calling.

I'm continually grateful to the community for all the contributions and
support. The [[mutect_commit][MuTect commit history for bcbio]] is a great example of multiple
distributed collaborators working towards the shared goal of integrating and
validating a caller.

#+LINK: bcbio http://github.com/chapmanb/bcbio-nextgen
#+LINK: dream https://www.synapse.org/#!Synapse:syn312572
#+LINK: dream_synthetic https://www.synapse.org/#!Synapse:syn312572/wiki/62018
#+LINK: dream_about https://bcbio-nextgen.readthedocs.org/en/latest/contents/testing.html#cancer-tumor-normal
#+LINK: jointval http://bcb.io/2014/10/07/joint-calling/
#+LINK: giab http://genomeinabottle.org/
#+LINK: icgc https://icgc.org/
#+LINK: icgc_pipelines http://biorxiv.org/content/early/2014/12/24/013177
#+LINK: icgc_calling http://biorxiv.org/content/early/2014/12/24/013177
#+LINK: bina http://www.bina.com/
#+LINK: bina_varsim https://github.com/bioinform/varsim
#+LINK: bina_varsim_cancer http://info.bina.com/hs-fs/hub/419826/file-1900487108-pdf/Posters/ASHG_2014_VarSim.pdf
#+LINK: bina_cancer http://info.bina.com/cancer-resources
#+LINK: justin https://twitter.com/BioInfo
#+LINK: miika https://github.com/mjafin
#+LINK: luca https://github.com/lbeltrame
#+LINK: az http://www.astrazeneca.com/Home
#+LINK: lorena https://github.com/lpantano
#+LINK: mutect_commit https://github.com/chapmanb/bcbio-nextgen/commits/master/bcbio/variation/mutect.py
#+LINK: wwcrc http://www.gla.ac.uk/researchinstitutes/cancersciences/ics/
#+LINK: hbc http://hsphbio.ghost.io/

* Variant caller validation

We used a large collection of open source variant callers to detect SNPs, Indels
and structural variants:

- [[mutect][MuTect (1.1.5)]] -- A SNP only caller, built on GATK UnifiedGenotyper, from the
  Broad Institute. MuTect requires a license if used for commercial purposes.
- [[vardict][VarDict (2014-12-15)]] -- A SNP and indel caller from Zhongwu Lai and the
  oncology team at AstraZeneca.
- [[freebayes][FreeBayes (0.9.20-1)]] -- A haplotype aware realigning caller for SNPs and
  indels from Erik Garrison and Gabor Marth's lab.
- [[varscan][VarScan (2.3.7)]] -- A heuristic/statistic based somatic SNP and indel caller
  from Dan Koboldt and The Genome Institute at Washington University.
- [[scalpel][Scalpel (0.3.1)]] -- Micro-assembly based Indel caller from Giuseppe Narzisi and
  the Schatz lab. We pair Scalpel with MuTect to provide a complete set of small
  variant calls.
- [[https://github.com/arq5x/lumpy-sv][LUMPY (0.2.7)]] -- A [[http://genomebiology.com/2014/15/6/R84/abstract][probabilistic structural variant caller]] incorporating both split
  read and read pair discordance, developed by Ryan Layer in
  [[http://quinlanlab.org/][Aaron Quinlan]] and [[http://genome.wustl.edu/people/groups/detail/hall-lab/][Ira Hall's]] labs.
- [[https://github.com/tobiasrausch/delly][DELLY (0.6.1)]] -- An [[http://bioinformatics.oxfordjournals.org/content/28/18/i333.abstract][integrated paired-end and split-end structural variant caller]]
  developed by Tobias Rausch.
- [[wham][WHAM (1.5.1)]] -- A structural variant caller that can incorporate association
  testing, from Zev Kronenberg in Mark Yandell's lab at the University of Utah.

bcbio runs these callers and uses simple ensemble methods to combine
[[ensemble_snp][small variants (SNPs, indels]]) and [[ensemble_sv][structural variants]] into final combined
callsets. The [[ensemble_bvr][new small variant ensemble method]] uses a simplified majority rule
classifier that picks variants to pass based on being present in
a configurable number of samples. This performs well and is faster than the
previous implementation that made use of both this approach and
a subsequent support vector machine step.

#+LINK: mutect https://www.broadinstitute.org/cancer/cga/mutect
#+LINK: vardict https://github.com/AstraZeneca-NGS/VarDict
#+LINK: varscan http://varscan.sourceforge.net/
#+LINK: scalpel http://scalpel.sourceforge.net/
#+LINK: freebayes https://github.com/ekg/freebayes
#+LINK: wham https://github.com/jewmanchue/wham
#+LINK: ensemble_snp http://bcb.io/2013/10/21/updated-comparison-of-variant-detection-methods-ensemble-freebayes-and-minimal-bam-preparation-pipelines/
#+LINK: ensemble_sv http://bcb.io/2014/08/12/validated-whole-genome-structural-variation-detection-using-multiple-callers/
#+LINK: ensemble_bvr https://github.com/chapmanb/bcbio.variation.recall

Using the 100x whole genome tumor/normal pair from [[dream_synthetic][DREAM synthetic dataset 3]]
we evaluated each of the callers for sensitivity and precision on small variants
(SNPs and indels). This synthetic dataset contains 100% tumor cellularity with 3
subclones at allele frequencies of 50%, 33% and 20%.

#+BEGIN_HTML
<a href="http://i.imgur.com/Ck8AJQk.png">
  <img src="http://i.imgur.com/Ck8AJQk.png" width="700"
       alt="Cancer evaluation for SNPs and indels">
</a>
#+END_HTML

In addition to the whole genome results, the [[cancer_val_album][validation album]] includes
results from running against the same dataset limited to exome regions. This has
identical patterns of sensitivity and precision. It runs quicker, so is useful
for evaluating changes to filtering or program parameters.

We also looked at structural variant calls for larger deletions, duplications and
inversions. Here is the precision and sensitivity for duplications across multiple
size classes:

#+BEGIN_HTML
<a href="http://i.imgur.com/csQPl5p.png">
  <img src="http://i.imgur.com/csQPl5p.png" width="700"
       alt="Cancer evaluation for structural variants -- duplications.">
</a>
#+END_HTML

The [[cancer_val_album][full album of validation results]] includes the comparisons for deletion and
inversion events. These comparisons measure the contributions of individual
callers within [[ensemblesv][an ensemble approach]] that attempts to maximize sensitivity and
specificity for the combined callset. Keep in mind that each of the individual
results make use of other caller information in filtering. Our goal is to create
the best possible combined calls, rather than a platform for unbiased
comparison of structural variant callers. We're also actively working on
improving sensitivity and specificity for individual callers and expect the
numbers to continue to evolve. For example, [[https://github.com/jewmanchue/][Zev Kronenberg]] added WHAM's ability
to identify different classes of structural changes, and we're still in the
process of improving the classifier.

#+LINK: cancer_val_album http://imgur.com/a/qba5k
#+LINK: ensemblesv http://bcb.io/2014/08/12/validated-whole-genome-structural-variation-detection-using-multiple-callers/

* Improvements in filtering

Our evaluation comparisons show best effort attempts to provide good quality
calls for every caller. The final results often come from multiple rounds of
improving sensitivity and precision by adjusting program parameters or
downstream filtering. The goal of tightly integrating bcbio with validation is
that the community can work on defining a set of parameters and filters that
work best in multiple cases, and then use these directly within the same
framework for processing real data.

In presenting the final results only, it may not be clear that plugging
a specific tool into a custom bash script will not always produce the same
results we see here. As an example, here are the improvements in FreeBayes
sensitivity and precision from our initial implementation, presented over the
exome regions of synthetic dataset 3:

#+BEGIN_HTML
<a href="http://i.imgur.com/NJFaoas.png">
  <img src="http://i.imgur.com/NJFaoas.png" width="500"
       alt="FreeBayes caller improvements">
</a>
#+END_HTML

The original implementation used a [[vcflib][vcfsamplediff based]] approach to filtering, as
recommended on the [[fb_filterorig][FreeBayes mailing list]]. The current, improved, version uses a
[[fb_filter][custom filter based on genotype likelihoods]], based on the approach in the
[[speedseq][speeseq pipeline]].

In general, you can find all of the integration work for individual callers
in the [[https://github.com/chapmanb/bcbio-nextgen/tree/master/bcbio/variation][bcbio source code]], broken down by caller. For instance, here is the
integration work on [[https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/variation/mutect.py][MuTect]]. The goal is to make all of the configuration
settings and filters fully transparent so users can understand how they work
when using bcbio, or transplant them into their own custom pipelines.

#+LINK: fb_filter https://github.com/chapmanb/bcbio-nextgen/blob/4fe770cc1343f8e1a3f3fab1771bad13eb94df7a/bcbio/variation/freebayes.py#L217
#+LINK: fb_filterorig https://groups.google.com/d/msg/freebayes/beLYRuHMkQE/RwFMniDmBYoJ
#+LINK: vcflib https://github.com/ekg/vcflib
#+LINK: speedseq https://github.com/cc2qe/speedseq
#+LINK: bcbiovdir

* Tumor-only prioritization

The above validations were all done on cancer calling with tumor and normal
pairs. The filters to separate pre-existing [[germline][germline mutations]] from
[[somatic][cancer specific somatic mutations]] rely on the presence of a normal
sample. In some cases, we don't have matched normal samples to do this filtering. Two
common examples are [[ffpe][FFPE samples]] and tumor cell lines. For these
samples, we'd like to be able to prioritize likely tumor specific variations for
followup using publicly available resources.

We implemented a prioritization strategy from tumor-only samples in bcbio that
takes advantage of publicly available resources like [[cosmic][COSMIC]], [[clinvar][ClinVar]],
[[1000g][1000 genomes]], [[esp][ESP]] and [[exac][ExAC]]. It uses [[gemini][GEMINI]] to annotate the initial tumor-only VCF
calls with external annotations, then extracts these to [[bcbio_priority][prioritize variants]] with
high or medium predicted impact, not present in 1000 genomes or ExAC at more
than 1% in any subpopulation, or identified as pathenogenic in COSMIC or ClinVar.

Validating this prioritization strategy requires real tumor samples with known
mutations. Our synthetic datasets are not useful here, since the variants
[[dream_synthetic][do not necessarily model standard biological variability]]. You could spike in
biologically relevant mutations, as done in the [[bina_varsim_cancer][VarSim cancer simulated data]],
but this will bias towards our prioritization approach since both would use the
same set of necessarily imperfect known variants and population level mutations.

We took the approach of using published tumor data with validated
mutations. [[https://twitter.com/druvus][Andreas Sjödin]] identified a [[val_paper][Hepatoblastoma exome sequencing paper]]
with [[val_dataset][publicly available sample data]] and 23 validated cancer related variations
across 5 samples. This is a baseline to help determine how stringent to be
in removing potential germline variants.

The prioritization enriches variants of interest by 35-50x without
losing sensitivity to confirmed variants:

| sample | caller    | confirmed | enrichment | additional | filtered |
|--------+-----------+-----------+------------+------------+----------|
| HB2T   | freebayes | 6 / 7     |        44x |       1288 |    56046 |
| HB2T   | mutect    | 6 / 7     |        48x |       1014 |    47755 |
| HB2T   | vardict   | 6 / 7     |        36x |       1464 |    52090 |
| HB3T   | freebayes | 4 / 4     |        46x |       1218 |    54997 |
| HB3T   | mutect    | 4 / 4     |        49x |        961 |    46894 |
| HB3T   | vardict   | 4 / 4     |        35x |       1511 |    51404 |
| HB6T   | freebayes | 4 / 4     |        43x |       1314 |    56240 |
| HB6T   | mutect    | 4 / 4     |        51x |        946 |    47747 |
| HB6T   | vardict   | 3 / 4     |        35x |       1497 |    51625 |
| HB8T   | freebayes | 6 / 6     |        42x |       1364 |    57121 |
| HB8T   | mutect    | 6 / 6     |        47x |       1053 |    48639 |
| HB8T   | vardict   | 6 / 6     |        35x |       1542 |    52642 |
| HB9T   | freebayes | 2 / 2     |        41x |       1420 |    57582 |
| HB9T   | mutect    | 2 / 2     |        44x |       1142 |    49858 |
| HB9T   | vardict   | 2 / 2     |        36x |       1488 |    53098 |

We consistently missed one confirmed mutation in the HB2T sample. This variant,
reported as a somatic mutation in an [[C2orf57][uncharacterized open reading frame (C2orf57]]),
may actually be a germline mutation in the study sub-population. The variant is
present at a [[missing_exac][10% frequency in the East Asian population]] but only 2% in the
overall population, based on data from both the ExAC and 1000 genomes
projects. Although the ethnicity of the original samples is not reported, the
study authors are all from China. This helps demonstrate the effectiveness of
large population frequencies, stratified by population, in prioritizing and
evaluating variant calls.

The major challenge with tumor-only prioritization approaches is that you can't
expect to accurately filter private germline mutations that you won't find in
genome-wide catalogs. With a sample set using a small number of validated
variants it's hard to estimate the proportion of 'additional' variants in
the table above that are germline false positives versus the proportion that are
additional tumor-only mutations not specifically evaluated in the study. We plan
to continue to refine filtering with additional validated datasets to help
improve this discrimination.

Practically, bcbio automatically runs prioritization with all tumor-only
analyses. It filters variants by adding a ``LowPriority`` filter to the output VCF
so users can readily identify variants flagged during the prioritization.

#+LINK: gemini https://github.com/arq5x/gemini
#+LINK: 1000g http://www.1000genomes.org/
#+LINK: esp http://evs.gs.washington.edu/EVS/
#+LINK: exac http://exac.broadinstitute.org/
#+LINK: clinvar http://www.clinvar.com/
#+LINK: cosmic http://cancer.sanger.ac.uk/cancergenome/projects/cosmic/
#+LINK: germline https://en.wikipedia.org/wiki/Germline_mutation
#+LINK: somatic https://en.wikipedia.org/wiki/Mutation#Somatic_mutations
#+LINK: ffpe http://www.mobio.com/blog/2012/07/20/formalin-fixed-paraffin-embedded-tissue-dna-isolation-the-basics/
#+LINK: bcbio_priority https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/variation/prioritize.py
#+LINK: dream_synthetic https://www.synapse.org/#!Synapse:syn312572/wiki/62018
#+LINK: val_paper http://onlinelibrary.wiley.com/doi/10.1002/hep.27243/full
#+LINK: val_dataset http://www.ebi.ac.uk/ena/data/view/SRP037747
#+LINK: missing_exac http://exac.broadinstitute.org/variant/2-232458818-C-T
#+LINK: C2orf57 http://useast.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000177673;r=2:231592901-231594283;t=ENST00000313965;redirect=no

* Future work

This is a baseline for assessing SNP, indel and structural variant calls in
cancer analyses. It also prioritizes impact variants in cases where we lack
normal matched normals. We plan to continue to improve cancer variant calling in
[[bcbio][bcbio]] and some areas of future focus include:

- Informing variant calling using estimates of tumor purity and sub-clonal
  frequency. bcbio integrates [[cnvkit][CNVkit]], a copy number caller, which exports
  read count segmentation data. Tools like [[theta][THetA2]], [[https://github.com/morrislab/phylowgs][phyloWGS]], [[https://github.com/uci-cbcl/PyLOH][PyLOH]]
  and [[https://github.com/genome/sciclone][sciClone]] use these inputs to estimate normal contamination and sub-clonal
  frequencies.

- Focusing on difficult to call genomic regions and provide additional
  mechanisms to better resolve and improve caller accuracy in these
  regions. Using small variant calls to define problematic genome areas with
  large structural rearrangements can help prioritize and target the use of more
  computationally expensive methods for copy number assessment, structural
  variant calling and de-novo assembly.

- Evaluating calling and tumor-only prioritization on
  [[horizon][Horizon reference standards]]. They contain a larger set of
  [[horizonvars][validated mutations]] at a variety of frequencies.

As always, we welcome suggestions, contributions and feedback.

#+LINK: horizon http://www.horizondx.com/our-products/q-seq-hdx-reference-standards
#+LINK: horizonvars http://www.horizondx.com/media/item/202
#+LINK: cnvkit  http://cnvkit.readthedocs.org/en/latest/
#+LINK: theta https://github.com/raphael-group/THetA
#+LINK: discovar http://www.broadinstitute.org/software/discovar/blog/
