An efficient algorithm for generating the internal branches of a Kingman coalescent

Abstract

Coalescent simulations are a widely used approach for simulating sample genealogies, but can become computationally burdensome in large samples. Methods exist to analytically calculate a sample's expected frequency spectrum without simulating full genealogies. However, statistics that rely on the distribution of the length of internal coalescent branches, such as the probability that two mutations of equal size arose on the same genealogical branch, have previously required full coalescent simulations to estimate. Here, we present a sampling method capable of efficiently generating limited portions of sample genealogies using a series of analytic equations that give probabilities for the number, start, and end of internal branches conditional on the number of final samples they subtend. These equations are independent of the coalescent waiting times and need only be calculated a single time, lending themselves to efficient computation. We compare our method with full coalescent simulations to show the resulting distribution of branch lengths and summary statistics are equivalent, but that for many conditions our method is at least 10 times faster.

Keywords: Coalescent; Coalescent simulations; Genealogical topology.

Figure 1

Density functions for the distribution of genealogical branch lengths with sizes between 2 and 10 for a sample of size 50

Figure 2

The exact probabilities of observing a given number of branches of a given size in the genealogy of a sample with size 50.

Figure 3

Topology free sampling. After calculating and storing intermediate topological probabilities, it is possible to efficiently generate segments of genealogies, and relevant summary statistics, using Monte Carlo integration.

Figure 4

(A) Cumulative distribution of branch lengths for branches of size 3, 5, and 7 in a genealogy from a sample with size 1,000. The curves give the probability of observing a branch shorter in length than the value given along the x-axis, with each color representing a different branch size. The double-dashed lines represents lengths drawn from 50,000 full coalescent simulations, the dashed lines from 50,000 genealogies sampled according to the topology free approach presented in the methods section, and the solid line from 50,000 genealogies sampled using expected rather than random waiting times. The results of the three approaches are so similar they are indistinguishable from the plots. (B) The number of branches with a size of 3, 5, or 7 in a genealogy from a sample with size 1,000. The solid black lines give the analytical values from formula (10) and the red bars were observed values from 50,000 genealogies realized with full coalescent simulations.

Figure 5

The sampling of individual branches allows us to calculate summary statistics with our topology free approach that are not available from the summed total length. (A) For sample sizes between 750 and 100 the variance between branch lengths was very similar between our topology free approach and full coalescent simulations. (B) In smaller sample sizes, larger size branches had higher inter-branch variance in coalescent simulations, due to the negative correlation between their lengths. (C) In larger samples the probability that two mutations of the same size arose on the same branch is nearly identical between coalescent simulations and topology free sampling. (D) In smaller sample sizes, the use of topology free sampling resulted in very similar, yet slightly lower, probability estimates than full coalescent simulations. This reflects the same underlying causes as seen in panel B, where our assumption that the number of branches is independent of their lengths and the allowing of shared coalescent events leading to lower inter-branch variance.