Fast simulation of identity-by-descent segments

Abstract

The worst-case runtime complexity to simulate haplotype segments identical by descent (IBD) is quadratic in sample size. We propose two main techniques to reduce the compute time, both of which are motivated by coalescent and recombination processes. We provide mathematical results that explain why our algorithm should outperform a naive implementation with high probability. In our experiments, we observe average compute times to simulate detectable IBD segments around a locus that scale approximately linearly in sample size and take a couple of seconds for sample sizes that are less than ten thousand diploid individuals. In contrast, we find that existing methods to simulate IBD segments take minutes to hours for sample sizes exceeding a few thousand diploid individuals. When using IBD segments to study recent positive selection around a locus, our efficient simulation algorithm makes feasible statistical inferences, e.g., parametric bootstrapping in analyses of large biobanks, that would be otherwise intractable.

Keywords: 60–08; 92D15; 92–04; 92–08; 92–10; coalescent; computational runtime; identity-by-descent.

Fig. 1

Conceptual framework for IBD segment lengths. (Left) Sample haplotypes $a$, $b$, $c$, $d$ trace their lineages back to common ancestors at times $t_4$, $t_4+t_3$, $t_4+t_3+t_2$. (Right) Relative to a focal point, the haplotype segments lengths $R_a$, $R_b$, $L_a$, $L_b$ are independent, identically distributed $\text{Exponential}(t_4)$. The lengths shared IBD are $R_{a,b}:=\text{min}(R_a,R_b)$ and $L_{a,b}:=\text{min}(L_a,L_b)$. The IBD segment length $W_{a,b}:=L_{a,b}+R_{a,b}\sim \text{Gamma}(2,2·t_4)$ exceeds the detection threshold $w$ Morgans.

Fig. 2

Compute time to simulate IBD segment lengths around a locus depending on algorithm implementation. Compute time (y-axis) in seconds by sample size (x-axis) in thousands is averaged over five simulations. The legend denotes colored line styles for implementations using Algorithm 1 as is (blue), merging only (orange), pruning only (green), and neither pruning nor merging (red). The main text describes “merging” and “pruning” techniques. The demography is the population bottleneck. The Morgans length threshold is 0.01.

Fig. 3

Compute time to simulate IBD segment lengths around a locus depending on the detection threshold and population size. Compute time (y-axis) in seconds by sample size (x-axis) in thousands is averaged over five simulations. The legends denote colored line styles for A) different detection thresholds (in Morgans) with $N={10}^5$ fixed or B) different population sizes with 0.02 Morgans fixed.

Fig. 4

Percentage of regression model predictions explained by linear and quadratic effects. The percentage of predicted compute time (y-axis) in seconds by sample size (x-axis) in thousands with respect to linear and quadratic effects. Plots show results for A) the constant population size $N\equiv N_e=\mathrm{10,\; 000}$ versus B) $N\equiv N_e=\mathrm{100,\; 000}$. The detectable IBD segments are simulated with a 0.02 Morgans threshold.