Sequencing a population at low coverage

A few reads per person cannot fix genotypes, and hard calls from them inflate rare variants and hide heterozygotes.

Read first: Molecular population genetics

Reads, not genotypes

A sequencer does not report genotypes. It reports reads, short pieces of DNA, and each read comes from one of the two chromosomes of one person. At a site, a person with mean coverage C = 4 has a number of reads that varies around 4 (the module uses a Poisson distribution). Each read picks one of the two chromosomes with probability ½, and each base in it is misread with probability ε, here 0.002.

Low depth hides heterozygotes. With d reads at a heterozygous site, the chance that every read comes from the same chromosome is 2 × (1/2)^d, which is 50% at d = 2 and 12.5% at d = 4. Those people look homozygous. A naive caller also needs a minimum depth, here 2 reads, so genotypes covered by 0 or 1 reads are missing. Calls then give too few heterozygotes and a lower heterozygosity.

Depth per genotype

Figure 1. Reads per genotype for 8 people sequenced at 4× over 3000 bp, seed 1. The ticks are the expectation. 9.8% of genotypes have fewer than 2 reads and are missing from the calls, against 9.2% expected.

Errors become singletons, likelihoods absorb them

The naive caller calls a heterozygote when the minority allele makes up at least τ = 0.2 of the reads. One misread base among 4 reads is 25%, so it becomes a false heterozygote. The same is true up to 5 reads. A sample of 8 people at 4× over 3000 bp expects 8 × 3000 × 4 × 0.002 = 192 misread bases. Each one that lands on a site with at most 5 reads creates a variant that is not there, and each such variant is seen once, a singleton.

Genotype likelihoods keep the uncertainty instead of calling. With 3 reference reads and 1 alternative read, the likelihood of the heterozygote is 94 times that of the reference homozygote. A single site cannot tell whether that is a true variant. The module calls a site variable only when the likelihood ratio against no variation reaches 23.9, which one stray read never does. It then estimates the spectrum without calling genotypes at all, by maximum likelihood (the GL (EM) estimate).

The module offers a middle option, GL calls, which test each site first and then call genotypes from the posterior. It removes the false variants, 0 in the sample below, but the test also rejects real variants that are seen in one or two reads. It finds 83 of the 100 true sites and 7 of the 19 singletons. Testing sites first trades false positives for missing rare variants, and the call-free estimate avoids that trade.

The spectrum from hard calls and from likelihoods

Figure 2. The site frequency spectrum of the same sample, projected to 12 chromosomes, seed 1. The sample has 100 true variable sites and 19 singletons. The naive calls find 204 variable sites and 108 singletons, and 108 sites with no true variant are called variable. The GL (EM) estimate has 98 sites and 18 singletons. Tajima's D is +0.92 in truth, −0.90 from naive calls, and +0.97 from GL (EM).

Bias in segregating sites against coverage

Figure 3. Estimated Watterson's θ divided by the true value, for 3 simulated samples of 8 people, 2000 bp, ε = 0.002. Each coverage re-sequences the same three samples, so only the reads change. At 4× the naive calls give 2.21 times the truth and the posterior GL calls give 0.71. At 16× all three methods give 1.00. GL (EM) stays between 0.91 and 1.02 at every coverage.

Bias in heterozygosity against coverage

Figure 4. Heterozygosity divided by the true value with ε = 0, so only missing and under-called heterozygotes matter. The naive calls give 0.57 at 1×, 0.80 at 4× and 0.98 at 16×. The dashed line is the closed-form expectation for the naive caller, 0.79 at 4×.

Exercise

Reproduce the default sample

Open the module, leave every setting at its default and press Run (the first run uses seed 20261). Read the θ_W and Tajima's D cards, and the false SNP count under the confusion matrix.

What to expect. The truth has 565 SNPs, θ_W = 140 and D = −0.33. The naive calls give θ_W = 315 and D = −1.73 and report 698 false SNPs. GL (EM) gives θ_W = 139 and D = −0.36. Hard calls make a constant-size sample look like a growing population.

Open the walkthrough step →

Exercise

Remove the errors, keep the low depth

Set Cov to 2 and Error to 0, then run. Compare the heterozygosity and θ_π cards.

What to expect. Naive heterozygosity falls to 0.0078 against a truth of 0.0112, about 30% too low. Naive θ_π barely moves, 128.8 against 128.4. Under-calling removes heterozygotes but leaves pairwise differences about right. At Cov = 30 all methods agree with the truth.

Open the walkthrough step →

What to keep

Reads are samples of chromosomes: a heterozygote covered by d reads looks homozygous with probability 2 × (1/2)^d.
Sequencing errors become false singletons in hard calls, which inflate θ_W and push Tajima's D negative, as growth would.
The naive θ_W bias peaks at moderate coverage (2.21 times at 4×) and vanishes at 16×, while naive heterozygosity stays under the truth at every coverage tested (0.98 at 16×).
The call-free GL (EM) estimate stays within 10% of the true θ_W and within 3% of the true heterozygosity at every coverage in these figures.