The effects of sampling strategy on the quality of reconstruction of viral population dynamics using Bayesian skyline family coalescent methods: A simulation study

Abstract

The ongoing large-scale increase in the total amount of genetic data for viruses and other pathogens has led to a situation in which it is often not possible to include every available sequence in a phylogenetic analysis and expect the procedure to complete in reasonable computational time. This raises questions about how a set of sequences should be selected for analysis, particularly if the data are used to infer more than just the phylogenetic tree itself. The design of sampling strategies for molecular epidemiology has been a neglected field of research. This article describes a large-scale simulation exercise that was undertaken to select an appropriate strategy when using the GMRF skygrid, one of the Bayesian skyline family of coalescent methods, in order to reconstruct past population dynamics. The simulated scenarios were intended to represent sampling for the population of an endemic virus across multiple geographical locations. Large phylogenies were simulated under a coalescent or structured coalescent model and sequences simulated from these trees; the resulting datasets were then downsampled for analyses according to a variety of schemes. Variation in results between different replicates of the same scheme was not insignificant, and as a result, we recommend that where possible analyses are repeated with different datasets in order to establish that elements of a reconstruction are not simply the result of the particular set of samples selected. We show that an individual stochastic choice of sequences can introduce spurious behaviour in the median line of the skygrid plot and that even marginal likelihood estimation can suggest complicated dynamics that were not in fact present. We recommend that the median line should not be used to infer historical events on its own. Sampling sequences with uniform probability with respect to both time and spatial location (deme) never performed worse than sampling with probability proportional to the effective population size at that time and in that location and frequently was superior. As a result, we recommend this approach in the design of future studies. We also confirm that the inclusion of many recent sequences from a single geographical location in an analysis tends to result in a spurious bottleneck effect in the reconstruction and caution against interpreting this as genuine.

Keywords: coalescent; phylodynamics; sampling; simulation.

Figure 1.

Depiction of the population structure used in structured coalescent simulations. Circles represent demes; two are small, two medium, and two large. Thick arrows represent fast rates of movement between demes (0.05 transitions per lineage per year) and thin arrows slower rates (0.025 per lineage per year).

Figure 2.

Overlaid median lines of reconstructed skygrid plots for 50 replicates of the uniform temporal sampling scheme in Scenario 1. The black line is the true population size.

Figure 3.

KDEs for the distribution of statistics indicating the accuracy and precision of the skygrid reconstructions in Scenario 1.

Figure 4.

Skygrid reconstructions for the 50 replicates of the uniform sampling scheme in Scenario 1, sorted by increasing percent error. The black line is the true EPS, the dark blue line the median estimate, and the 95 per cent HPD region is in light blue.

Figure 5.

Scatter plots depicting the relationship between sample size and statistics used to evaluate the performance of the skygrid reconstruction for 100 replicates of the uniform sampling scheme in Scenario 1. The red line represents the best-fit model determined by weighted least squares regression and corrected Akaike information criterion and the grey lines the limit of the 95 per cent confidence interval given by the fitted values of the error term of this model. (a) plots adjusted percent error against sample size and (b) HPD size against sample size.

Figure 6.

KDEs for the distribution of statistics indicating the accuracy and precision of the skygrid reconstructions in Scenario 2. Distributions are coloured by sampling scheme.

Figure 7.

KDEs for the distribution of statistics indicating the accuracy and precision of the skygrid reconstructions in Scenario 3. Distributions are coloured by sampling scheme.

Figure 8.

Overlaid median lines of 50 reconstructed skygrid plots for Scenario 5, where additional samples are selected from one deme in the last 0.25 years of the timeline. The black line is the true population size. Plot titles refer to the oversampled deme and its size.

Figure 9.

KDEs for the distribution of statistics indicating the accuracy and precision of the skygrid reconstructions in Scenario 4. Each graph plots the distributions of each statistic for several temporal sampling schemes when the spatial sampling scheme is fixed.

Figure 10.

KDEs for the distribution of statistics indicating the accuracy and precision of the skygrid reconstructions in Scenario 4. Each graph plots the distributions of each statistic for several spatial sampling schemes when the temporal sampling scheme is fixed.

Figure 11.

KDEs for the distribution of statistics indicating the accuracy and precision of the skygrid reconstructions in Scenario 5. Distributions are coloured by sampling scheme.