A family-based probabilistic method for capturing de novo mutations from high-throughput short-read sequencing data

Abstract

Recent advances in high-throughput DNA sequencing technologies and associated statistical analyses have enabled in-depth analysis of whole-genome sequences. As this technology is applied to a growing number of individual human genomes, entire families are now being sequenced. Information contained within the pedigree of a sequenced family can be leveraged when inferring the donors' genotypes. The presence of a de novo mutation within the pedigree is indicated by a violation of Mendelian inheritance laws. Here, we present a method for probabilistically inferring genotypes across a pedigree using high-throughput sequencing data and producing the posterior probability of de novo mutation at each genomic site examined. This framework can be used to disentangle the effects of germline and somatic mutational processes and to simultaneously estimate the effect of sequencing error and the initial genetic variation in the population from which the founders of the pedigree arise. This approach is examined in detail through simulations and areas for method improvement are noted. By applying this method to data from members of a well-defined nuclear family with accurate pedigree information, the stage is set to make the most direct estimates of the human mutation rate to date.

Figure 1

Under-calling heterozygous genotypes affects de novo detection at a given site. In the top panel, the mother’s genotype is called AA by the sample-independent approach since the binomial probability of sampling once the C allele among 20 reads if the mother is heterozygote is very small. (Nielsen et al., 2011, suggest calling a site homozygous if the minor allele is less than 20%, a rule which we adopt for these examples.) When the family data is considered jointly, identifying a C in the child increases the probability of the AC genotype for the mother, leading to a low probability of de novo mutation at this site. (It is much more likely that the mother’s chromosomes were sampled unevenly, ≈10⁻⁵, than that there is an actual mutation at the site, ≈10⁻⁸.) In the bottom panel, the child’s genotype is called AA. However, given an error rate and the parental coverage, the probability of a de novo mutation at this site is high. The de novo mutation probabilities were computed using the method described here with the following parameters: θ = 0.001, ε = 0.005, μ = μ_s = 2 × 10⁻⁸. (See section 2.1 for a description of these parameters.)

Figure 2

Trio model for a single site. Nucleotide bases from the sequencing reads aligning to the site of interest (R = {R_M, R_F, R_O}) are the observed data. Neither individual genotypes nor their transmission pattern are observed; they are the hidden data. Double lines denote the transmission/sampling of diploid genotypes, while single lines haploid genotypes. The parental zygotic genotypes (m = m_am_b and f = f_a f_b) are sampled from a common population and are the founding alleles of the pedigree. Wavy lines denote where sequencing takes place. The parameters in the model have been placed in proximity to the branches that they affect.

Figure 3

Transition points in model inferences. (A) transition from error to de novo mutation, (B) effect of depth on error-mutation transition, (C) effect of low parent depth, (D) effect of low offspring depth, (E) transition from de novo mutation to inherited allele, (F) transition from de novo mutation to inferred inheritance. The parameter values used are θ = 0.001, ε = 0.01, μ = 2 × 10⁻⁷, μ_s = 0.0. The read structures are given in the tables; each row represents the read data taken for a family member, with the offspring data (O) as the bottom row, and the columns represent the number of reads observed for each allele. ‘G’ and ‘T’ read counts are always 0.

Figure 3

Transition points in model inferences. (A) transition from error to de novo mutation, (B) effect of depth on error-mutation transition, (C) effect of low parent depth, (D) effect of low offspring depth, (E) transition from de novo mutation to inherited allele, (F) transition from de novo mutation to inferred inheritance. The parameter values used are θ = 0.001, ε = 0.01, μ = 2 × 10⁻⁷, μ_s = 0.0. The read structures are given in the tables; each row represents the read data taken for a family member, with the offspring data (O) as the bottom row, and the columns represent the number of reads observed for each allele. ‘G’ and ‘T’ read counts are always 0.

Figure 4

Rate estimates from repeated trio simulations. Ten simulation replicates of trios with 20× coverage were generated. Parameters were estimated for each replicate. ‘Simulated’ rates were calculated from the full data, and ‘estimated’ rates from the observable data. Simulated rates only differ from simulation parameters due to the stochastic nature of the simulations. Both sets are plotted relative to the simulation parameters: μ = 10⁻⁶, ε = 0.01007, and θ = 0.001.

Figure 5

Trio Simulation ROC curves. For each simulated set, the corresponding ROC curve and AUC value are presented based on all calls with δ ≥ 0.01. The parameter values used in simulations are θ = 0.001, ε = 0.01, μ = 1 × 10⁻⁶, and μ_s = 0.0. The dashed line shown along the diagonal represents the expected ROC curve for a random, un-useful classifier. A perfect classifier goes straight up then straight across. The “-Ctrl” results are for simulations in which all reads are perfectly aligned back to the reference genome.

Figure 6

De novo mutations predictions from trio simulations. 259 sites contained de novo mutations in the simulations (gray line). The total number of mutation calls at three different levels of δ are given by the triangles for each simulation. The circles indicate the amount of true positives for each δ threshold. The distances between triangles and circles represent the amount of false positives, and the distances between the circles and gray line gives the number of false negatives. The “-Ctrl” columns contain results for simulations in which all reads are perfectly aligned back to the reference genome.

Figure 7

Validation results stratified by de novo calling algorithm. Three methods were used to produce candidate sites of de novo mutation in two families (Conrad et al., 2011, Durbin et al., 2010). These candidates were then experimentally validated and classified as either “true mutations” (germline, somatic, or cell-line), “false mutations” (inherited genetic variation or no variation at all), or “inconclusive”. Sites not reported by a method but reported by another method appear in the No Call bar. FPIR corresponds to the method described in this paper. Modified from Conrad et al. (2011).

Figure A.1

Possible genealogical trees for a sample of four alleles from a single population. The times of coalescent events, t_i, are expessed in 2N_e generations, where N_e is the effective population size. Type I occurs twice as often as Type II.