BCFtools/trio-dnm
Petr Danecek
Version: November 25, 2020
This document describes methods implemented in the bcftools/trio-dnm plugin. At the time of writing,the plugin implements two de novo mutation (DNM) discovery models: the original ”genotype-centric”
approach, which is also at the core of the Family-aware Illumina Genotype Likelihood-based method(FIGL) and is used by DeNovoGear (DNG) [Conrad et al., 2011, Ramu et al., 2013], and a new ”allelecentric” approach which aims to overcome some shortcomings and bugs in DeNovoGear implementation,and which is in these notes referred to as ”trio-dnm model”.

1

Genotype-centric approach (DeNovoGear model)

The sole purpose of reimplementing the DeNovoGear model was to build trust in the new method byshowing that BCFtools can successfully reproduce results from the original model. This method can be
invoked by running the plugin with the --use DNG option.
The DNG model relies on genotype likelihoods contained in the FORMAT/PL annotation calculated bybcftools/mpileup (originally samtools/mpileup). It is beyond the scope of this document to explain thisstep in full detail, therefore we only show the final formula and refer the interested reader to [Li, 2010]and htslib/errmod.c code.
For simplicity of notation we will assume bi-allelic sites, with x, y ∈ {0, 1} representing the reference(0) and the alternate (1) allele, and consider only diploid genotypes, i.e. G ∈ {xx, xy, yy}. Assume weobserved in total n reads covering a genomic position, of which k were alternate reads and n − k reference
reads. Each of the bases is assigned the error probability i derived from the base quality and mapping
quality.
Genotype likelihoods express how consistent is the observed data D with all possible genotypes
P (D|00) =k−1
Yfi
βni(i ),

i=0

P (D|11) =n−1Yfi
βni(i ),i=k nP (D|01) =1/2nk(because   1/2),

fi
where βni
is the conditional probability of observing more than i errors while accounting for their possible
dependency using the empirical parameter fi = 0.97η l + 0.03.

These are the input data for the DNG model, for more details see the section ”Revised MAQ model”
in [Li, 2010]. In addition DeNovoGear requires FORMAT/DP annotation to filter low coverage sites(DP < 10), but since it is not used in the calculation, it is not required by bcftools/trio-dnm.
The DNG model evaluates joint data likelihoods for all possible combinations of trio genotypes P , M ,F of the proband, mother and father as follows
LP,M,F = P (D|P, M, F ) · P (P |M, F ) · P (M, F ).
1(1)The first term represents the product of genotype likelihoods of the proband, mother and father; the second term represents the transmission probability (0.25, 0.5, or 1 if compatible with Mendelian inheritanceor µ = 10−8 for each mutated allele); and the third term is the prior probability of drawing two genotypes
M and F from the population. For details see the supplementary materials of [Conrad et al., 2011] andbcftools/trio-dnm code.
The most likely combination of trio genotypes P, M, F incompatible with Mendelian inheritance is thenselected and the posterior probability calculated as
LP,M,F,
P (is DNM) = X
Lp,m,fp,m,fwhere summation in the denominator is over all possible genotype combinations, including combinations compatible with Mendelian inheritance. The program adds the rounded phred-scaled value of
1 − P (is DNM) as the FORMAT/DNM annotation, bigger values indicate increased likelihood of thevariant being a true de novo mutation.
In a test performed on 30,898 candidate de novo mutations there was perfect correlation between resultsfrom DeNovoGear and bcftools/trio-dnm implementations (R2 = 1.0). Note that the DNG implementation in bcftools/trio-dnm uses the same mutation rates for both SNPs and indels and does not attempt
to reproduce the calculation for sex chromosomes for reasons given further.

Figure 1. DNM probability as calculated by the original DeNovoGear program and its
reimplementation in BCFtools/trio-dnm.

2Allele-centric approach (trio-dnm model)

The DeNovoGear model suffers from a shortcoming which manifests itself in inflated false positive rate.
The underlying reason is that genotype likelihoods were designed for germline variant calling and, to adegree, are robust to mapping and alignment errors. Specifically, if the observed fraction of alternatebases deviates strongly from the expected binomial sampling probability, the likelihood of heterozygous
genotype will be smaller than the likelihood of homozygous genotype despite strong presence of thealternate allele. This is a property desired for germline variant calling but is unsuitable for DNM
discovery. Indeed, the most common failure mode of the DNG model is the situation where alternate
reads are present in both the proband and a parent, but as most likely are deemed the heterozygous
genotype in the proband and homozygous genotypes in both parents. Although this behavior may leadto serendipitous discovery of valid biological cases (true mosaic mutations in parents), such discoveries
are accidental and indistinguishable from false calls on the basis of call quality.

2The revised approach described in this section fixes this problem by sensing the presence of alternate
alleles in the parents directly, rather than relying on genotype likelihoods. Genotype likelihoods are used
only to determine heterozygous genotype in the proband.

2.1

Allele Quality

Let x,i denote the base error as before (i.e. the maximum of base and mapping error), but this timefor each base x ∈ {A, C, G, T } separately. For simplicity of notation we limit the description to single
nucleotides, but indels are handled the same way. We introduce a new variable ax which will represent
the probability of base x being genuine and not an artifact. Then if the base x is present in reads coveringa genomic position kx times, the probability ax can be expressed as(Qkxx,i if kx > 0,1 − i=0(2)ax =0otherwise.
This quantity is calculated at the bcftools/mpileup step by adding the -a QS option and it will be used
to calculate a modified trio genotype likelihood by replacing P (D|P, M, F ) in equation (1) with a new00term P 0 (D|Pxy , Mxy,uv, Fxy,rs) as explained next.
Extend the notation of genotype likelihood in the proband by adding subscripts to indicate the genotype
in question, and of the parental likelihoods by adding the allelic context of the proband’s genotype. For
0example, M01,00
stands for the likelihood of the homozygous reference genotype in the mother (0/0)
conditioned on the heterozygous genotype in the proband (0/1). The modified trio genotype likelihood
can be then expressed as0000
P 0 (D|Pxy , Mxy,uv, Fxy,rs) = Pxy · Mxy,uv· Fxy,rswith maternal genotype likelihoods calculated as4 
miY0Mxy,uv=1 − mi
i=0 1 − miif i ∈ {u, v},if i ∈ {u, v} \ {x, y},otherwise, if u = v,(3)(4)

where mx denotes the allele quality ax in mother (and similarly for father). The first case checks fora strong presence of the alleles in the parental genotype uv; the second penalizes the presence of theassumed de novo allele in the parental genotype; and the third penalizes presence of multiple alleles inhomozygous parental genotypes.

2.2

Modifications to DeNovoGear calculations

The math notes in [Conrad et al., 2011] and the DeNovoGear code contain some errors and in some cases
lack clarity.
1. The math notes state that conditional on a segregating site with 1 or 2 alternate alleles, genotypes
with 1 alternate allele are found with frequency 3/5 and genotypes with 2 alternate alleles are foundwith frequency 2/5. This is in fact reversed, there are 4 configurations with a single alternate allele
and 6 with two alternate alleles, giving frequencies 2/5 and 3/5, respectively:
00 0100 1001 00
10 00.. 4 out of 10 = 2/500 1101 01
10 01
01 1010 1011 00.. 6 out of 10 = 3/53

2. It is not clear how the term 0.001 (as opposed to 0.002) has been derived in the math notes.
3. The DeNovoGear code deviates from the math notes, setting for example the prior probability ofa single alternate allele in parents as
g_priors[i] = 0.995 * 0.002 * (3.0 / 5.0) * (4.0 / 5.0) * 0.5;

or triallelic cases as
g_priors[i] = 0.002 * 0.002 / 414;

No explanation is provided as for where the terms 0.995, 0.5 or 414 come from.
4. The DNG code for setting chrX priors has multiple bugs. For example, the transmission of chrXfrom mother to a boy does not depend on the paternal genotype, etc. From this reason bcftools/triodnm does not attempt to reimplement this part of the code.

2.3

Revised DeNovoGear priors

Label the reference allele 0 and consider at most three non-reference alleles x ∈ {1, 2, 3}. Set the priorprobabilities of sampling reference and alternate genotypes as follows
P1 = 0.998
the probability of sampling two reference genotypes at sites with one alternate allele common
in the population (biallelic sites)
P3 = (1 − P1 )2
the probability of sampling a non-reference genotype at sites with two different alternate alleles
common in the population (triallelic sites)
P 2 = 1 − P1 − P3 ≈ 1 − P1
the probability of sampling a non-reference genotype at biallelic sites
P4 = 10−26
the probability of observing tetra-allelic site
2.3.1

Autosomal chromosomes

The prior probability of sampling four chromosomes from the population for parental genotypes M andF is then as follows
0 alternate alleles:P1M, F = 00, 00




1 alternate allele:
P2 · (1/15) · (1/3)
M, F ∈ {xx, xx} where x ∈ {1, 2, 3}P2 · (2/15) · (1/3)
M, F ∈ {00, xx ; xx, 00}P2 · (4/15) · (1/3)
M, F ∈ {0x, xx ; x0, xx ; xx, 0x ; xx, x0}
P (M, F ) =P2 · (4/15) · (1/3)
M, F ∈ {0x, 0x ; 0x, x0 ; x0, 0x ; x0, x0}





P2 · (4/15) · (1/3)
M, F ∈ {00, 0x ; 00, x0 ; 0x, 00 ; x0, 00}






2 alternate alleles:P3 · (1/19) · (1/3)
M, F ∈ {00, xy ; 0x, xy ; . . . }

3 or more alternate alleles:


P44

2.3.2Chromosome X in males

For chromosome X inheritance pattern in males we define analogously
PX,1 = 0.999
PX,3 = (1 − PX,1 )2
PX,2 = 1 − PX,1 − PX,3 ≈ 1 − PX,1The priors for drawing the maternal genotype then become
0 alternate alleles:PX,1M = 00
1 alternate allele:
PX (M ) =PX,2 · (1/3) · (1/3)
M = xx where x ∈ {1, 2, 3}P·(2/3)·(1/3)M ∈ {0x ; x0}
X,2





2 alternate alleles:PX,3 · (1/3)
M ∈ {xy ; xz ; yz}
2.3.3

Chromosome X in females

For chromosome X inheritance in females we approximate as
0 alternate alleles:P1M, F = 00, 0





1 alternate allele:
P2 · (3/7) · (1/3)
M, F ∈ {00, x ; 0x, 0 ; x0, 0} where x ∈ {1, 2, 3}
PXX (M, F ) =P·(3/7)·(1/3)
M, F ∈ {0x, x ; 0x, x ; xx, 0}
2



P2 · (1/7) · (1/3)
M, F = xx, x
2 alternate alleles:P3 · (1/9) · (1/3)
M, F ∈ {0x, y ; . . . ; xx, y ; . . . ; xy, y ; . . . }5

3Results

Overall performance
The test data consisted of 50,764 candidate calls from de novos.ddd 10k.29-01-2018.tab.gz. For theseplots calls on chrX were included, we note that the graphs look the same with chrX excluded.

ts/tv (mark at 15e3 DNMs)

2.65

1e4
3.0

BCFtools/trio-dnm
DeNovoGear

Number of DNMs (mark at tstv=2.4)

2.70
2.60
2.55
2.50
2.45
2.40
2.35
2.30

BCFtools/trio-dnm
DeNovoGear

2.5
2.0
1.5
1.0
0.5

20

40

60

80
Threshold

100

120

140

2040

60

80
Threshold100

120140Figure 2. The ts/tv ratio and the number of sites at increasing DNM quality score cutoff. Circles
mark points in the graph with the same number of sites and ts/tv, respectively. The qualities
calculated by the new model (BCFtools/trio-dnm) are placed higher than those calculated by
DeNovoGear and thus contain less noise in the same-sized call set.

1e4
2.6

BCFtools/trio-dnm
DeNovoGear

4.5

2.5
Number of DNMs

4.0

2.4
ts/tv

BCFtools/trio-dnm
DeNovoGear

2.3
2.2

3.5
3.0
2.5

2.1

2.0

2.0

1.5
0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
Sensitivity

0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
Sensitivity

Figure 3. The ts/tv ratio and the number of sites at increasing DNM quality score cutoff plotted as afunction of sensitivity with respect to a clean set of 15,213 high-quality de novos from coding/splicing
regions. Note that the truth set was generated by filtering the DeNovoGear calls and therefore it
should have higher sensitivity.

6

Cases missed by DeNovoGear
In all IGV screenshots that follow, displayed are from top to bottom father, mother, and the proband.
The blue and pink color shows the forward and reverse strand orientation of the reads.

Figure 4. Call missed by DeNovoGear, chromosome X, the proband is male.

7

Figure 5. Call missed by DeNovoGear, chromosome X, the proband is female.

Figure 6. Call missed by DeNovoGear, chromosome 1.
False calls by DeNovoGear

Figure 7. This is likely a false call caused by mapping artifact, the alternate allele appears in bothparents. Can’t exclude parental mosaicity.
8

Cases missed by BCFtools/trio-dnm

Figure 8. Arguably, this might be a valid biological case with low-level mosaicity in the father.
However, it could be also argued that the program correctly determined that the alternate allele is not
de novo in the proband and a different calling model should be used. Note that BCFtools/trio-dnm has
a parameter that allows to adjust noise tolerance and effectively call cases like this at the cost ofdecreased specificity.

Uncertain cases - should these be called or not?

Figure 9. Mosaic variant in the child or an artifact?

9

Figure 10. Is this a valid DNM? Note the single alternate read in the father.

10References[Conrad et al., 2011] Conrad D, et al. Variation in genome-wide mutation rates within and between
human families Nat Genet, 43, 712-714 (2011).
[Ramu et al., 2013] Ramu A, et al. DeNovoGear: de novo indel and point mutation discovery and phasing, Nat Methods, 10(10), 985-987 (2013).
[Li, 2010] Li H, The revised MAQ model, (2010)
http://samtools.github.io/bcftools/samtools.pdf

11