Joint genotype inference with germline and somatic mutations

The joint sequencing of related genomes has become an important means to discover rare variants. Normal-tumor 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 cancer-related death among children. With the aim of uncovering the full spectrum of germline and somatic genetic alterations in childhood ALL genomes, we conducted whole-exome re-sequencing on a unique cohort of over 120 exomes of childhood ALL quartets, each comprising a patient's tumor and matched-normal 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 normal-tumor 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 open-source 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 cancer-related death among children. Advances in the understanding of the pathobiology of ALL have led to risk-targeted 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 so-called 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.
Next-generation sequencing (NGS) technologies [1] have enabled the genome-wide identification of human disease-related variants. Analysis pipelines have been established for the large-scale 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 normal-tumor pairs.
In general, genetic relationships (like normal-tumor 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 normal-tumor 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 target-enrichment techniques [8] are available for sequencing only regions of interest. In particular, sequencing the socalled exome covering all gene-coding regions, has been a routine step in medical applications [9]. Through our ongoing paediatric oncogenomics study, we conducted whole exome deep re-sequencing of a unique cohort of over 120 exomes of childhood ALL quartets, consisting of the patient's tumor and matched-normal 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 normal-tumor samples, including SomaticSniper [14], MutationSeq [15], and JointSNVMix [16]. SomaticSniper and MutationSeq employ machine-learning 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 tumor-normal 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 Figure 1 Alignments for the normal, tumor, and parental reads samples highlight putative somatic and germline mutations. The illustration shows a region along chromosome 12, displayed in the Integrative Genomics Viewer [7]. Mismatched bases are highlighted. that genotypes and lineage-specific 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 tumor-sampling impurity and sequencing errors. All model parameters are estimated in an expectation-maximization algorithm [18]. Figure 2 illustrates our model of single-nucleotide polymorphisms at homologous loci across four genomes linked by inheritance and somatic mutations. The model quantifies the descent-by-modification 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.

Probabilistic model
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.  computed by using a standard DNA substitution model quantifying the divergence from the reference genome assembly. Assuming the simple Jukes-Cantor model with a reference nucleotide genotype y, π [x] = ν/3 for x = y and Π y = 1 − v, where v is the parental genome's divergence from the reference. More general divergence models and known SNP frequencies can be accommodated by using π [x] = y∈{A,C,G,T} π ref y v y→x, where V x→y is the model's substitution probability from x to y and ref is a known allele frequency.

Parental genotype priors
For a diploid locus, let φ (xx) and φ xy denote the frequency of homozygous xx and heterozygous xy genotypes. Allele frequencies immediately determine diploid allele frequencies in standard Hardy-Weinberg equilibrium: φ(xx) = (π [x]) 2 for a homozygous locus xx, and φ(xy) = 2π [x]π [y] for a heterozygous xy . In order to accommodate different homozygous-heterozygous ratios, we employ an additional parameter γ ≤ 1.
(In other words, g is numerically analogous to the socalled inbreeding coefficient, or F-statistic.) 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 (g F , g M ), 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 v x→x = P X = x|X = x . Let χ (g F , g M → g N ) denote the probability of normal genotype g N given the mutated parental genotypes g F , g M . 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 basecall-error probability pairs (z k , Î k ) for k = 1, ..., m. The k-th 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 BAM-format file gives the Sanger-encoded sequencing error in the QUAL column or the OQ tag.) Namely, where q k is the reported quality score, and φ = 10 1/10 = 0.794 · · ·. 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 ρ[y] = P{y k = y} for each possible allele y. At a homozygous locus xx, r[y] = 1 and r[ Impure tumor samples have a mixed distribution, which is the linear combination of the normal and tumor genotype distributions. For tumor reads, r[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 basecall-error 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 r is defined as Since base calls are independent across reads when conditioned on the allele mixture in the sample,

Complete likelihood
Let g F, g M, g N, g Tdenote the diploid genotypes for father, mother, normal, and tumor samples, respectively. Let g F and g M 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 R consists of base call and error pairs (z k , Î k ). The likelihood for the aligned base calls is then The L(g F , ...) 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 r that are determined by genotypes and tumor sampling purity, as discussed above (see Allele sampling and sample impurity).
The p(g F , ...) 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 Expectation-Maximization (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 χ (g F , g M → g N ) as 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 b = r[z] appears in the formula. Define If no reads call z, then T b [z] = 1. Equation (10)  Impure tumor sample. The pure tumor (g T ) and normal genotypes (g N ) 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 r, 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 g N → g T. There are 6 possible queried allele frequencies b ≠ 0, 1, 1/2 ( Figure 3): 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 sample-specific H' 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 R at a given locus for a given sample (by allele and quality score), and use run-length encoding to achieve a compact characterization h(R). The encoding is not used for higher coverages or widely varying quality scores, where h would take too many bits. Information-theoretic considerations [21] suggest that compactly encoded R occur more often ( h is our proxy for Kolmogorov complexity). Short codes h are placed in a small hash table to find recurrent calls (in our experiments with 20-30× total coverage by AB SOLiD sequencing reads, 20-30% 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 prob- Model parameters are optimized using the EM algorithm [18]. In one iteration, likelihoods and various posterior probabilities are computed across all loci in the so-called E-step, which are then used to set the model parameters for the next iteration in the so-called M-step. 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 Jukes-Cantor model for the germline and tumor mutations is fairly straightforward by using posterior probabilities. For example, the tumor mutation parameter ν (T) is set by summing the posterior probabilities for allele substitutions across all loci: where N is the number of loci, p j (g N , g T ) is the posterior probability for a normal-tumor genotype pair at locus j, and a(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 b values) by reading the input read-mapping files at each iteration. At the same time, we compute a calibrated map μ : {0, 1, . . . , 93} → [0, 1] from reported base-calling qualities to sequencing errors in the same framework. Note that both μ and ω have well-estimated initial values (μ starts with the canonical Phred-scaled 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 g /2 probability each. Otherwise, with probability (1 -g), the output is the same heterozygote xy. Homozygous genotypes are output without any change. Clearly, if the input genotype distribution is for Hardy-Weinberg, then the machine's output is distributed by the probabilities of (1). Accordingly, the divergence and heterozygosity parameters for the parental genotype prior j 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 g.

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 Sainte-Justine 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.

Whole-genome sequencing
Tumor and normal DNA samples from Quartet B were submitted to Illumina, Inc, for deep whole-genome sequencing using the standard operating procedures of the HiSeq 2000 sequencing platform. Table 1 summarizes coverage statistics and tumor impurity. The whole-genome 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 parental-offspring 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 sample-specific 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 E-step of the optimization, where sums of posterior probabilities are calculated across loci, is well-suited 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.

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 5-6× 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 (wall-clock time) on standard multi-core computer workstations with 16 Gbytes of memory.
Exome-sequencing data from Quartet A was used to assess the concordance of genotyping calls by QUADGT and a well-established variant caller, The Genome Analysis Toolkit [12]. We used independently produced whole-genome (WG) sequencing reads for the normaltumor pair in Quartet B (with 124× total coverage, see Table 1) to gauge the two variant callers' sensitivity. 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.

Sensitivity assessment
The normal-tumor pair in Quartet B was submitted to Illumina, Inc. for deep whole-genome (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 non-reference 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 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 Expectation-Maximization) is close to the Illumina's estimate from whole-genome data (see Table 1). Table 4 suggests that the joint variant-calling 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.    Table 4 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 normal-tumor 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, transition-transversion ratios and nucleotide composition, as well as site-specific priors relying on public variant databases and gene annotations.  Somatic mutation calls of QUADGT (qGT columns) and the Genome Analysis Toolkit (GATK columns) were ranked by quality score (smaller of normal and tumor genotype call qualities in GATK's VCF output, and score of somatic call in QUADGT's VCF output that is q = -10 log 10 (1 -p) for a mutation with posterior probability p). Rank is "-" at loci where no somatic mutations were inferred by QUADGT or GATK. Genome coordinates refer to NCBI 36.1/hg18 genome build.