Abstract
The joint sequencing of related genomes has become an important means to discover rare variants. Normaltumor genome pairs are routinely sequenced together to find somatic mutations and their associations with different cancers. Parental and sibling genomes reveal de novo germline mutations and inheritance patterns related to Mendelian diseases.
Acute lymphoblastic leukemia (ALL) is the most common paediatric cancer and the leading cause of cancerrelated death among children. With the aim of uncovering the full spectrum of germline and somatic genetic alterations in childhood ALL genomes, we conducted wholeexome resequencing on a unique cohort of over 120 exomes of childhood ALL quartets, each comprising a patient's tumor and matchednormal material, and DNA from both parents. We developed a general probabilistic model for such quartet sequencing reads mapped to the reference human genome. The model is used to infer joint genotypes at homologous loci across a normaltumor genome pair and two parental genomes.
We describe the algorithms and data structures for genotype inference, model parameter training. We implemented the methods in an opensource software package (QUADGT) that uses the standard file formats of the 1000 Genomes Project. Our method's utility is illustrated on quartets from the ALL cohort.
Background
Acute lymphoblastic leukemia (ALL) is the most common paediatric cancer and the leading cause of cancerrelated death among children. Advances in the understanding of the pathobiology of ALL have led to risktargeted treatment regimes and increased survival rates, but treatment is still far from optimal. Childhood ALL arises after the acquisition of a series of DNA sequence abnormalities. These initiating events, or socalled driver mutations, ultimately confer a selective growth advantage, and are causally implicated in cancer development. A central goal of cancer genome analysis is the identification of cancer genes that, by definition, carry driver mutations.
Nextgeneration sequencing (NGS) technologies [1] have enabled the genomewide identification of human diseaserelated variants. Analysis pipelines have been established for the largescale sequencing of individual tumor genomes [2]. Briefly, short sequencing reads are collected from the tumor sample, mapped to the reference genome assembly, and the set of aligned reads are used to infer variations across the genome at homologous loci covered with multiple reads. The sequence variants in the tumor genome may be the result of somatic mutations, or constitutional variants preserved in the somatic lineage. In order to distinguish somatic mutations from conserved variants, it is necessary to sequence normal and tumor samples side by side. The Cancer Genome Atlas Network [3] catalogues somatic mutations in different cancers using such normaltumor pairs.
In general, genetic relationships (like normaltumor pairs) can be efficiently exploited in genotype inference [4,5]. Inherited and de novo mutations can be traced through jointly sequenced family relatives [6]. Here, we consider variant detection in normaltumor pairs coupled with parental samples. Such quartet data are used to categorize variants in the tumor and normal genomes by their origin: see Figure 1. One can readily classify inherited variants and de novo germline mutations by comparing the genotypes in the trio of normal and parental genomes. Likewise, somatic mutations correspond to differences between the normal and tumor genomes.
While sequencing an entire human genome is still too expensive for the average research laboratory, various targetenrichment techniques [8] are available for sequencing only regions of interest. In particular, sequencing the socalled exome covering all genecoding regions, has been a routine step in medical applications [9]. Through our ongoing paediatric oncogenomics study, we conducted whole exome deep resequencing of a unique cohort of over 120 exomes of childhood ALL quartets, consisting of the patient's tumor and matchednormal material as well as DNA from both parents.
Existing software tools
Various bioinformatics tools have been developed for genotyping individual genomes from sequencing data, including SNVMix [10], VarScan [11], and The Genome Analysis Toolkit GATK [12,13]. A couple of methods have been developed for the purpose of joint genotyping of paired normaltumor samples, including SomaticSniper [14], MutationSeq [15], and JointSNVMix [16]. SomaticSniper and MutationSeq employ machinelearning techniques for variant classification; JointSNVMix is based on a full Bayesian model incorporating prior genotype distributions, somatic mutations, and sequencing base call errors. The Strelka software package [17] infers joint tumornormal genotypes in a Bayesian model that also considers tumor sampling impurity: DNA collected from the tumor sample is usually "contaminated" to some degree with the normal tissue, and therefore the sequencing reads come from a mixture of normal and tumor genomes. To our knowledge, no existing variant caller incorporates somatic and germline mutation models simultaneously to handle quartet data as in our data sets.
Our contribution
We infer the four genotypes jointly in a framework that respects the rules of inheritance in the germline and somatic lineages. Aside from assigning belief to de novo and somatic mutations, we hypothesized, constrained patterns in one lineage have an indirect beneficial effect on the inference in other lineages. In particular, the "triangulation" of the normal genome by related genomes means that genotypes and lineagespecific mutations can be resolved more reliably: information from the parental genotypes reinforce the inference of somatic mutations, and tumor sequencing reads help to recognize constitutional mutations. We present a Bayesian framework that incorporates prior parental genotypes, inherited, de novo and somatic mutations, as well as tumorsampling impurity and sequencing errors. All model parameters are estimated in an expectationmaximization algorithm [18].
Methods
Probabilistic model
Figure 2 illustrates our model of singlenucleotide polymorphisms at homologous loci across four genomes linked by inheritance and somatic mutations. The model quantifies the descentbymodification relationships between the unknown genotypes via three sets of parameters. First, a genotype frequency model is assumed for the parental genotypes. Second, we assume a standard DNA substitution model for the frequency of germline mutations. The parental diploid genotypes determine the child's normal genotype by Mendelian inheritance. (For simplicity, we discuss only diploid genotypes: our implementation considers sex chromosomes in an analogous manner, but using the appropriate inheritance model.) Finally, another DNA substitution model with its own parameters determines mutations in the tumor genome.
Figure 2. Probabilistic model for base calls covering a single locus, with dependencies among genotypes and sequencing reads. See Equation (5) for the corresponding likelihood formula.
Base calls are assumed to be independent between different loci. The input for genotype inference at a single locus consists of nucleotide base calls made with their accompanying sequencing error probabilities.
Parental genotype priors
Let π[x] denote the frequency of each allele x ∈ {A; C; G; T} at a given locus. The prior allele frequencies are computed by using a standard DNA substitution model quantifying the divergence from the reference genome assembly. Assuming the simple JukesCantor model with a reference nucleotide genotype for and , where v is the parental genome's divergence from the reference. More general divergence models and known SNP frequencies can be accommodated by using , where is the model's substitution probability from to and is a known allele frequency.
For a diploid locus, let denote the frequency of homozygous and heterozygous genotypes. Allele frequencies immediately determine diploid allele frequencies in standard HardyWeinberg equilibrium: for a homozygous locus , and for a heterozygous . In order to accommodate different homozygousheterozygous ratios, we employ an additional parameter .
(In other words, γ is numerically analogous to the socalled inbreeding coefficient, or Fstatistic.) It is easy to verify that haploid allele frequencies are the same at any γ setting.
Mutations and inheritance
The child's normal genotype is determined by Mendelian inheritance and de novo point mutations with probabilities v_{x}_{→y }that occur within the parental germlines. For simplicity, germline mutations in a parental lineage (g_{F, }g_{M }for father and mother, respectively) are conceived of as mutations that result in a diploid genotype, which then determine the child's normal genotype by Mendel's laws.
Germline mutations follow standard molecular evolution model for substitutions in DNA. Let X denote the parent's normal allele at a locus, and X' denote the same allele at the end of the germline before gametogenesis. The mutation model specifies the probabilities that apply to every locus . Let denote the probability of normal genotype g_{N }given the mutated parental genotypes . Then χ may be 0, 1, 1/2 or 1/4, depending on the common alleles between the three genotypes.
Tumor genotype
The tumor genome undergoes mutations following the same type of molecular evolution model as the one used for germline mutations, but has its own parameters.
Sequencing errors
The alignment at the locus is represented as a set of basecallerror probability pairs (z_{k}, ∈_{k}) for k = 1, ..., m. The kth read calls allele z_{k }with an error probability of 0 <∈_{k }≤ 1. Typically, the aligned sequencing reads give nucleotide base calls and accompanying error probabilities on a logarithmic integer scale [19,20], called the Phred scale. (In principle, the input SAM or BAMformat file gives the Sangerencoded sequencing error in the QUAL column or the OQ tag.) Namely,
where q_{k }is the reported quality score, and . Let Z denote the base call, and Y denote the true nucleotide. We assume that errors are unbiased in the sense that
Allele sampling and sample impurity
Aligned sequencing reads randomly sample the haploid alleles at a given locus. Let y_{k }be the true allele for base call z_{k}. The locus' diploid genotype determines the frequency for each possible allele y. At a homozygous locus xx, ρ[y] = 1 and ρ[x'] = 0 for all x' ≠ x. At a heterozygous locus xx', ρ[y] = 1/2, ρ[x'] = 1/2 and ρ[x"] = 0 for all x" ≠ x, x'.
Impure tumor samples have a mixed distribution, which is the linear combination of the normal and tumor genotype distributions. For tumor reads, ρ[y] comes from a mixture (0 ≤ ω ≤ 1; ω = 1 for pure tumor sampling) between the normal and tumor genotypes:
Likelihood for aligned reads given the genotypes
Suppose we are given the set of basecallerror probability pairs (z_{k}, ∈_{k}) for k = 1, ..., m, representing the alignment at a locus. Let Z_{k }be the random variable for base call in read k at a fixed error rate ∈_{k}, and let Y_{k }be the random variable for the true sampled allele. Define the base call probability
The read likelihood for a given allele distribution ρ is defined as
Since base calls are independent across reads when conditioned on the allele mixture in the sample,
Complete likelihood
Let , , , denote the diploid genotypes for father, mother, normal, and tumor samples, respectively. Let and denote the parental genotypes after germline mutations. These six random variables constitute the hidden variables in our probabilistic model. The input is a set of aligned sequenced bases from each of the four samples:, , , . Every set consists of base call and error pairs (z_{k}, ∈_{k}). The likelihood for the aligned base calls is then
The factor is the likelihood for the reads, given the genotypes. By the independence of the sequencing runs,
The four factors are defined by (4), via the allele frequencies ρ that are determined by genotypes and tumor sampling purity, as discussed above (see Allele sampling and sample impurity).
The factor in (5) covers all mutation and inheritance events, as well as the parental genotypes. By the dependencies depicted in Figure 2,
Algorithmic techniques and data structures
Our algorithmic solutions address the efficient calculation of the likelihood formula of (5), and its use in an ExpectationMaximization (EM) framework for model parameter setting. First, we examine a straightforward decomposition of the likelihood formula dictated by the assumed probabilistic graphical model.
For the EM algorithm, we need to recompute likelihoods and posterior probabilities in a number of iterations, which can be directly achieved by storing all sequencing reads in memory, but such an approach may be costly. We scrutinize the computation of read likelihoods, in order to arrive at an economical data structure, also discussed in some detail, that eliminates the need to store all base calls in memory.
Likelihood decomposition
The summation formula for the full likelihood in Equation (5) is rearranged for efficiency, using the independencies apparent in (6) and (7). In addition, the germline mutations can be combined with Mendelian inheritance: define as
Then
If there are G possible diploid genotypes (G = 10 for DNA with four alleles), Equation (8) shows that the likelihood can be computed in O(G^{3}) time, instead of O(G^{6}) suggested by the definition of (5). In particular, the likelihood computation proceeds by calculating the following values.
Read likelihoods
Equation (4) is conveniently rearranged by different base calls:
Each subproduct is calculated separately as
by Equation (2). Note that only the called base's frequency β = ρ[z] appears in the formula. Define
If no reads call z, then T_{β}[z] = 1. Equation (10) becomes L_{z}(ρ) = T_{ρ}_{[z]}[z]. For pure diploid samples, ρ[z] may be 0, 1 or 1/2, corresponding to the possible subproducts for diploid samples E[y] = T_{0}[y] (sequencing error), C[y] = T_{1}[y] (homozygote yy), and H[y] = T_{1/2}[y] (heterozygote with y). Likelihood formulas become even more economical with normalized subproducts C'[y] = C[y]/E[y], H'[y] = H[y]/E[y], and scaling factor
Homozygous sample. For a pure homozygous sample with genotype yy (ρ[y] = 1 and ρ[z] = 0 for z ≠ y),
Heterozygous sample. For a pure heterozygous sample yy' (y ≠ y'; ρ[y] = ρ[y'] = 1/2 and ρ[z] = 0 for z ≠ y, y'), the likelihood becomes
Impure tumor sample. The pure tumor and normal genotypes are proper diploid genotypes. Tumor sequencing reads come from an impure sample: they sample the tumor genotype with probability ω, and the normal genotype with probability (1  ω). By Eq. (3), identical genotypes correspond to identical allele frequencies ρ, no matter what the purity level ω is. Suppose, however, that the locus has a mixture of divergent normal and tumor genotypes. Figure 3 shows that there are up to four correct base calls appearing in the reads, depending on the tumor mutation pattern . There are 6 possible queried allele frequencies β ≠ 0, 1, 1/2 (Figure 3):
Figure 3. Allele frequencies with impure normaltumor genotypes that differ from each other. Reads sample the tumor with probability ω, and the normal genotype with 1  ω.
If the normal and tumor genotypes are identical, then the purity is immaterial. The formulas with somatic mutations are:
Data structure for storing base calls
Storing individual base calls at each locus is costly, because the sequencing coverage may be large across the four samples. It suffices, however, to store the partial likelihood factors appearing in Equations (11), (12) and (14). In particular, for each locus with mapped base calls, it is enough to store the samplespecific H'[x], C'[x], and values for each called base x and the six possible values of β, in addition to a single scaling value E. For a sample with m base calls at the locus, (1 + 8m) variables thus suffice, independently of read coverage. The stored partial likelihoods are reused throughout the iterations optimizing the model parameters at a fixed tumor purity level ω.
Recurrent base calls
In our experience, loci with identical sets of base calls reoccur at an appreciable frequency, especially at lower coverages (less than about 15×) that characterize exon boundaries in exome sequencing. We exploit recurrent patterns in the piledup base calls to achieve even better memory usage and speed. Namely, we sort the base calls at a given locus for a given sample (by allele and quality score), and use runlength encoding to achieve a compact characterization . The encoding is not used for higher coverages or widely varying quality scores, where h would take too many bits. Informationtheoretic considerations [21] suggest that compactly encoded occur more often ( is our proxy for Kolmogorov complexity). Short codes are placed in a small hash table to find recurrent calls (in our experiments with 2030× total coverage by AB SOLiD sequencing reads, 2030% savings can be achieved this way).
Parameter optimization and genotype inference
The genotypes, germline and tumor mutations are inferred by carrying out the summation of (5) for a restricted set of genotypes in order to calculate posterior probabilities. For example, in order to infer the child's normal genotype, calculate first
Then the child's normal genotype has posterior probabilities .
Model parameters are optimized using the EM algorithm [18]. In one iteration, likelihoods and various posterior probabilities are computed across all loci in the socalled Estep, which are then used to set the model parameters for the next iteration in the socalled Mstep. The iterations continue until convergence is achieved. Among the optimized model parameters, the parental genotype priors, the germline and tumor mutation parameters are optimized through multiple iterations using the same set of precalculated partial likelihoods.
Setting the single parameter of the JukesCantor model for the germline and tumor mutations is fairly straightforward by using posterior probabilities. For example, the tumor mutation parameter is set by summing the posterior probabilities for allele substitutions across all loci:
where N is the number of loci, is the posterior probability for a normaltumor genotype pair at locus j, and α(g, g') is the expected number of substitutions for the two alleles given the diploid genotypes.
In order to set the tumor purity ω, the partial likelihoods need to be recomputed (for different T_{β }values) by reading the input readmapping files at each iteration. At the same time, we compute a calibrated map from reported basecalling qualities to sequencing errors in the same framework. Note that both μ and ω have wellestimated initial values (μ starts with the canonical Phredscaled values, and ω is estimated experimentally).
Decomposing zygosity and divergence
For the purposes of parameter inference, consider the following machine realizing the formulas for the parental genotype priors of (1). Upon receiving a heterozygous genotype xy, it flips either allele to output homozygous xx or yy with γ /2 probability each. Otherwise, with probability (1  γ), the output is the same heterozygote xy. Homozygous genotypes are output without any change. Clearly, if the input genotype distribution is for HardyWeinberg, then the machine's output is distributed by the probabilities of (1). Accordingly, the divergence and heterozygosity parameters for the parental genotype prior ϕ are inferred by treating the machine's input genotype as a hidden variable. Expected frequencies for divergent input genotypes are used to estimate the divergence parameter, and expected frequencies for heterozygous → homozygous "flips" are used to estimate γ.
Sequencing data
Exome sequencing
We conducted validation experiments using exomesequencing reads for two sets of quartets (A and B) generated on the Child Health Genomics Platform of the SainteJustine UHC Research Center. Sequencing reads were produced on Applied Biosystem's SOLiD sequencer and mapped with the accompanying software. Table 1 summarizes statistics on the mapped sequencing reads.
Table 1. Coverage statistics for two quartets used in validation experiments.
Wholegenome sequencing
Tumor and normal DNA samples from Quartet B were submitted to Illumina, Inc, for deep wholegenome sequencing using the standard operating procedures of the HiSeq 2000 sequencing platform. Table 1 summarizes coverage statistics and tumor impurity. The wholegenome data was further analyzed for somatic mutations with CASAVA and Strelka [17] by Illumina, Inc.
Implementation
We incorporated the presented methods into an opensource Java software package called QUADGT, using the standard file formats of the 1000 Genomes Project (SAM v1.4 [20] for input and VCF v4.1 [22] for output). Any of the input files may be missing, which makes QUADGT suitable to analyze sets with just normaltumor samples, or just parentaloffspring trios.
The probabilistic framework enabled us to couple the inference with confidence measures in the form of quality scores computed from posterior probabilities. The quality scores accompany samplespecific genotype calls (VCF's GQ field), as well as the joint genotyping for the four samples (VCF's QUAL column). The posterior probabilities for germline and somatic mutations are calculated, as well, by summing across all pertinent genotype assignments. The program specifically introduces ambiguity in the genotyping calls to meet prescribed quality scores: a definite x/y base call in a sample is replaced by x/. or ./y, and then by ./. in a greedy manner, in order to achieve high specificity.
Parallelization
The Estep of the optimization, where sums of posterior probabilities are calculated across loci, is wellsuited to parallelization. In our implementation, we use a hash key computed at each locus to assign the calculations to different computing threads running in parallel.
Availability
The QUADGT software package is publicly available at http://www.iro.umontreal.ca/~csuros/quadgt/ webcite.
Results and discussion
We used two quartet data sets (A and B) to compare independent and joint variant detection. Figure 1 summarizes the coverage statistics for the quartets. Exomesequencing reads at 56× coverage per sample were mapped to the human reference, and we inferred genome variants using our software package QUADGT. The entire analysis pipeline for one quartet set, including model training and genotyping, took about 12 hours (wallclock time) on standard multicore computer workstations with 16 Gbytes of memory.
Exomesequencing data from Quartet A was used to assess the concordance of genotyping calls by QUADGT and a wellestablished variant caller, The Genome Analysis Toolkit [12]. We used independently produced wholegenome (WG) sequencing reads for the normaltumor pair in Quartet B (with 124× total coverage, see Table 1) to gauge the two variant callers' sensitivity.
Concordance experiments
Table 2 compares individual genotyping calls made by the Genome Analysis Toolkit [12] and QUADGT on a small example consisting of calls on chromosome 12 (parameters were set to result in a comparable number of genotyping calls). The table illustrates that most calls are made in agreement between the two programs. The known relationships between the samples ensure the consistency of calls made by QUADGT, resulting in only 4 putative de novo germline mutations. GATK, ignorant of the relations, has 327 cases where a normal allele does not appear at either parent, which is by at least two magnitudes higher than what one would expect based on human intergeneration mutation rates [6]. GATK genotypes imply a large number (520) of somatic mutations, as well. As with de novo mutations, the joint calls by QUADGT are more conservative: only 14 somatic mutation calls are made.
Table 2. Comparison of calls made by the Genome Analysis Toolkit and QUADGT on Quartet A.
Sensitivity assessment
The normaltumor pair in Quartet B was submitted to Illumina, Inc. for deep wholegenome (WG) sequencing and somatic mutation calling. Table 3 tallies the WG somatic calls with coverage by exome data. Based on the WG sequencing data, 1817 loci have somatic mutations, of which 40 are covered by exome reads to sufficient depth (Table 3, b). Some of the 40 WG somatic calls have no or weak support (lines c and e) in exome reads, since with at most one exception, all normal and tumor base calls are identical with the reference. The remaining 24 WG somatic calls (line f) with sufficient exome read coverage and nonreference base calls are not beyond the reach of the variant callers to discover somatic mutations. On closer inspection, 5 WG somatic calls fall into a 150 bp region (line g), which, judged by largely divergent base calls, is likely to be misaligned; 4 WG somatic calls (line h) may even be erroneous. For instance, at chr8:10078796, which has a fairly low WG somatic quality score, the parental exome reads suggest that all four samples contain a heterozygous SNP that is in fact found in dbSNP (rs112078536). The remaining 15 WG somatic calls (line i) have support in the exome reads, and QUADGT discovers them all at some quality threshold cutoffs.
Table 3. Wholegenome somatic loci and exome genotyping on Quartet B.
Table 4 compares the sensitivity of joint and independent genotyping using the 15 WG somatic calls with support in exome base calls (Table 3, line i). First, it is notable that QUADGT's tumor purity estimation (by ExpectationMaximization) is close to the Illumina's estimate from wholegenome data (see Table 1).
Table 4. Exome somatic calls supported by wholegenome data.
Table 4 suggests that the joint variantcalling in QUADGT leads to better sensitivity than GATK, which does not consider the relations between the genomes. In particular, 9 out of 276 (3%) SOMATIC calls by QUADGT with quality score at least 30 have support in the WG data, whereas only 7 of GATK's 667 (1%) divergent genotypes of same quality threshold are validated. At a quality score cutoff of 20, 12 out of 555 (2%) and 7 out of 1312 (0.5%) are validated QUADGT and GATK predictions.
Conclusions
Sequencing multiple genomes with known pedigrees or clonal relationships has a great promise for understanding the development of particular diseases. Our experiments with sequenced quartets of parents and normaltumor pairs illustrate that the joint calls improve the reliability of inferred de novo and somatic mutations. The constraints imposed by the known relationships greatly improve the consistency of calls between different samples, and ultimately help to delineate the single nucleotide polymorphisms that can be associated with the disease.
A future release of the software is now under development that incorporates more nuanced substitution models with variable rates, transitiontransversion ratios and nucleotide composition, as well as sitespecific priors relying on public variant databases and gene annotations.
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
The authors are indebted to the patients and their parents for participating in this study. This study was supported by research funds provided by the Terry Fox Research Institute, the Canadian Institutes for Health Research, and Canada's National Sciences and Engineering Research Council. JFS is the recipient of a Cole Foundation studentship. RV is the recipient of a postdoctoral research fellowship from the Government of Canada through the "Foreign Affairs and International Trade Canada." DS holds the FrançoisKarlViau Research Chair in Pediatric Oncogenomics.
Declarations
Funding for the publication of this article was provided by a grant from the National Sciences and Engineering Research Council.
This article has been published as part of BMC Bioinformatics Volume 14 Supplement 5, 2013: Proceedings of the Third Annual RECOMB Satellite Workshop on Massively Parallel Sequencing (RECOMBseq 2013). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/S5.
References

Shendure J, Li H: Nextgeneration DNA sequencing.
Nature Biotechnology 2008, 26(10):11351145. PubMed Abstract  Publisher Full Text

Wood LD, Parsons DW, Jones S, Lin J, Sjöblum T, et al.: The genomic landscapes of human breast and colorectal cancers.
Science 2007, 318(5853):11081113. PubMed Abstract  Publisher Full Text

The Cancer Genome Atlas Research Network: Comprehensive genomic characterization defines human glioblastoma genes and core pathways.
Nature 2008, 455:10611068. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Le SQ, Durbin R: SNP detection and genotyping from lowcoverage sequencing data on multiple diploid samples.
Genome Research 2011, 21(6):952960. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Li B, Chen W, Zhan X, Busonero F, Sanna S, Sidore C, Cucca F, Kang HM, Abecasis GR: A likelihoodbased framework for variant calling and De Novo mutation detection in families.
PLoS Genetics 2012, 8(10):e1002944. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Roach JC, Glusman G, Smit AFA, Huff CD, Hubley R, Shannon PT, Rowen L, Pant KP, Goodman N, Bamshad M, Shendure J, Drmanac R, Jorde LB, Hood L, Galas DJ: Analysis of genetic inheritance in a family quartet by wholegenome sequencing.
Science 2010, 328(5978):636639. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP: Intergrative genomics viewer.
Nature Biotechnology 2011, 29:2426. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mamanova L, Coffey AJ, Scott CE, Kozarewa I, Turner EH, Kumar A, Howard E, Shendure J, Turner DJ: Targetenrichment strategies for nextgeneration sequencing.
Nature Methods 2010, 7:111118. PubMed Abstract  Publisher Full Text

Teer JK, Mullikin JC: Exome sequencing: the sweet spot before whole genomes.
Human Molecular Genetics 2010, 19:R145R151. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Goya R, Sun MG, Morin RD, Leung G, Ha G, Wiegand KC, Senz J, Crisan A, Marra MA, Hirst M, Huntsman D, Murphy KP, Aparicio S, Shah SP: SNVMix: predicting single nucleotide variants from nextgeneration sequencing of tumors.
Bioinformatics 2010, 26(6):730736. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Koboldt DC, Chen K, Wylie T, Larson DE, McLellan MD, Mardis ER, Weinstock GM, Wilson RK, Ding L: VarScan: variant detection in massively parallel sequencing of individual and pooled samples.
Bioinformatics 2009, 25(17):22832285. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA: The Genome Analysis Toolkit: A MapReduce framework for analyzing nextgeneration DNA sequencing data.
Genome Research 2010, 20(9):12971303. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

DePristo MA, et al.: A framework for variation discovery and genotyping using nextgeneration DNA sequencing data.
Nature Genetics 2011, 43:491498. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Larson DE, Harris CC, Chen K, Koboldt DC, Abbott TE, Dooling DJ, Ley TJ, Mardis ER, Wilson RK, Ding L: SomaticSniper: identification of somatic point mutations in whole genome sequencing data.
Bioinformatics 2012, 28(3):311317. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ding J, Bashashati A, Roth A, Oloumi A, Tse K, Zeng T, Haffari G, Hirst M, Marra MA, Condon A, Aparicio S, Shah SP: Feature based classifiers for somatic mutation detection in tumournormal paired sequencing data. [http:/ / bioinformatics.oxfordjournals.org/ content/ early/ 2011/ 11/ 13/ bioinformatics.btr629.abstract] webcite

Roth A, Ding J, Morin R, Crisan A, Ha G, Giuliany R, Bashashati A, Hirst M, Turashvili G, Oloumi A, Marra MA, Aparicio S, Shah SP: JointSNVMix: a probabilistic model for accurate detection of somatic mutations in normal/tumour paired nextgeneration sequencing data.
Bioinformatics 2012, 28(7):907913. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Saunders CT, Wong WSW, Swamy S, Becq J, Murray LJ, Cheetham RK: Strelka: accurate somatic smallvariant calling from sequenced tumornormal sample pairs.
Bioinformatics 2012, 28(14):18111817. PubMed Abstract  Publisher Full Text

Dempster AP, Laird NM, Rubin DP: Maximum likelihood from incomplete data via the EM algorithm.
Journal of the Royal Statistical Society Series B 1977, 39:138.

Ewing B, Green P: Basecalling of automated sequencer traces using Phred: II. error probabilities.
Genome Research 1998, 8:186194. PubMed Abstract  Publisher Full Text

Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Subgroup GPDP: The Sequence Alignment/Map format and SAMtools.
Bioinformatics 2009, 25(16):20782079. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Li M, Vitányi P: An Introduction to Kolmogorov Complexity and Its Applications. 3rd edition. Springer Science+Business Media; 2008.

Danecek P, Auton A, et al.: The variant call format and VCFTools.
Bioinformatics 2011, 27(15):21562158. PubMed Abstract  Publisher Full Text  PubMed Central Full Text