Evaluation of methods for the inference of ancestral recombination graphs

The ancestral recombination graph (ARG) is a structure that describes the joint genealogies of sampled DNA sequences along the genome. Recent computational methods have made impressive progress towards scalably estimating whole-genome genealogies. In addition to inferring the ARG, some of these methods can also provide ARGs sampled from a defined posterior distribution. Obtaining good samples of ARGs is crucial for quantifying statistical uncertainty and for estimating population genetic parameters such as effective population size, mutation rate, and allele age. Here, we use simulations to benchmark three popular ARG inference programs: ARGweaver, Relate, and tsdate. We use neutral coalescent simulations to 1) compare the true coalescence times to the inferred times at each locus; 2) compare the distribution of coalescence times across all loci to the expected exponential distribution; 3) evaluate whether the sampled coalescence times have the properties expected of a valid posterior distribution. We find that inferred coalescence times at each locus are more accurate in ARGweaver and Relate than in tsdate. However, all three methods tend to overestimate small coalescence times and underestimate large ones. Lastly, the posterior distribution of ARGweaver is closer to the expected posterior distribution than Relate's, but this higher accuracy comes at a substantial trade-off in scalability. The best choice of method will depend on the number and length of input sequences and on the goal of downstream analyses, and we provide guidelines for the best practices.

Bayesian Inference of Joint Coalescence Times for Sampled Sequences

The site frequency spectrum (SFS) is a commonly used statistic to summarize genetic variation in a sample of genomic sequences from a population. Such a genomic sample is associated with an imputed genealogical history with attributes such as branch lengths, coalescence times and the time to the most recent common ancestor (TMRCA) as well as topological and combinatorial properties. We present a Bayesian model for sampling from the joint posterior distribution of coalescence times conditional on the SFS associated with a sample of sequences in the absence of selection. In this model, the combinatorial properties of a genealogy, which is represented as a coalescent tree, are expressed as matrices. This facilitates the calculation of likelihoods and the effective sampling of the entire space of tree structures according to the Equal Rates Markov (or Yule-type) measure. Unlike previous methods, assumptions as to the type of stochastic process that generated the genealogical tree are not required. Novel approaches to defining both uninformative and informative prior distributions are employed. The uncertainty in inference due to the stochastic nature of mutation and the unknown tree structure is expressed by the shape of the posterior distributions. The method is implemented using the general purpose Markov Chain Monte Carlo software PyMC3. From the sampled posterior distribution of coalescence times, one can also infer related quantities such as the number of ancestors of a sample at a given time in the past (ancestral distribution) and the probability of specific relationships between branch lengths (for example, that the most recent branch is longer than all the others). The performance of the method is evaluated against simulated data and is also applied to historic mitochondrial data from the Nuu-Chah-Nulth people of North America. The method can be used to obtain estimates of the TMRCA of the sample. The relationship of these estimates to those given by ''Thomson's estimator'' is explored. Keywords: coalescent theory; Bayesian inference; time to most recent common ancestor; site frequency spectrum

Heterogeneity in effective size across the genome: effects on the Inverse Instantaneous Coalescence Rate (IICR) and implications for demographic inference under linked selection.

The relative contribution of selection and neutrality in shaping species genetic diversity is one of the most central and controversial questions in evolutionary theory. Genomic data provide growing evidence that linked selection, i.e. the modification of genetic diversity at neutral sites through linkage with selected sites, might be pervasive over the genome. Several studies proposed that linked selection could be modelled as first approximation by a local reduction (e.g. purifying selection, selective sweeps) or increase (e.g. balancing selection) of effective population size (Ne). At the genome-wide scale, this leads to a large variance of Ne from one region to another, reflecting the heterogeneity of selective constraints and recombination rates between regions. We investigate here the consequences of this variation of Ne on the genome-wide distribution of coalescence times. The underlying motivation concerns the impact of linked selection on demographic inference, because the distribution of coalescence times is at the heart of several important demographic inference approaches. Using the concept of Inverse Instantaneous Coalescence Rate, we demonstrate that in a panmictic population, linked selection always results in a spurious apparent decrease of Ne along time. Balancing selection has a particularly large effect, even when it concerns a very small part of the genome. We quantify the expected magnitude of the spurious decrease of Ne in humans and Drosophila melanogaster, based on Ne distributions inferred from real data in these species. We also find that the effect of linked selection can be significantly reduced by that of population structure.

Spectral analysis of transient amplifiers for death–birth updating constructed from regular graphs

A central question of evolutionary dynamics on graphs is whether or not a mutation introduced in a population of residents survives and eventually even spreads to the whole population, or becomes extinct. The outcome naturally depends on the fitness of the mutant and the rules by which mutants and residents may propagate on the network, but arguably the most determining factor is the network structure. Some structured networks are transient amplifiers. They increase for a certain fitness range the fixation probability of beneficial mutations as compared to a well-mixed population. We study a perturbation method for identifying transient amplifiers for death–birth updating. The method involves calculating the coalescence times of random walks on graphs and finding the vertex with the largest remeeting time. If the graph is perturbed by removing an edge from this vertex, there is a certain likelihood that the resulting perturbed graph is a transient amplifier. We test all pairwise nonisomorphic regular graphs up to a certain order and thus cover the whole structural range expressible by these graphs. For cubic and quartic regular graphs we find a sufficiently large number of transient amplifiers. For these networks we carry out a spectral analysis and show that the graphs from which transient amplifiers can be constructed share certain structural properties. Identifying spectral and structural properties may promote finding and designing such networks.

A geometrical framework for f –statistics

A detailed derivation of the f–statistics formalism is made from a geometrical framework. It is shown that the f–statistics appear when a genetic distance matrix is constrained to describe a four population phylogenetic tree. The choice of genetic metric is crucial and plays an outstanding role as regards the tree–like–ness criterion. The case of lack of treeness is interpreted in the formalism as presence of population admixture. In this respect, four formulas are given to estimate the admixture proportions. One of them is the so–called f4–ratio estimate and we show that a second one is related to a known result developed in terms of the fixation index FST. An illustrative numerical simulation of admixture proportion estimates is included. Relationships of the formalism with coalescence times and pairwise sequence differences are also provided.

Influence of the Emulsifier System on Breakup and Coalescence of Oil Droplets during Atomization of Oil-In-Water Emulsions

Spray drying of whey protein-based emulsions is a common task in food engineering. Lipophilic, low molecular weight emulsifiers including lecithin, citrem, and mono- and diglycerides, are commonly added to the formulations, as they are expected to improve the processing and shelf life stability of the products. During the atomization step of spray drying, the emulsions are subjected to high stresses, which can lead to breakup and subsequent coalescence of the oil droplets. The extent of these phenomena is expected to be greatly influenced by the emulsifiers in the system. The focus of this study was therefore set on the changes in the oil droplet size of whey protein-based emulsions during atomization, as affected by the addition of low molecular weight emulsifiers. Atomization experiments were performed with emulsions stabilized either with whey protein isolate (WPI), or with combinations of WPI and lecithin, WPI and citrem, and WPI and mono- and diglycerides. The addition of lecithin promoted oil droplet breakup during atomization and improved droplet stabilization against coalescence. The addition of citrem and of mono- and diglycerides did not affect oil droplet breakup, but greatly promoted coalescence of the oil droplets. In order to elucidate the underlying mechanisms, measurements of interfacial tensions and coalescence times in single droplets experiments were performed and correlated to the atomization experiments. The results on oil droplet breakup were in good accordance with the observed differences in the interfacial tension measurements. The results on oil droplet coalescence correlated only to a limited extent with the results of coalescence times of single droplet experiments.