Estimation of population growth or decline in genetically monitored populations

Effects of Social Structure on Effective Population Size Change Estimates

Most methods currently used to infer the "demographic history of species" interpret this expression as a history of population size changes. The detection, quantification, and dating of demographic changes often rely on the assumption that population structure can be neglected. However, most vertebrates are typically organized in populations subdivided into social groups that are usually ignored in the interpretation of genetic data. This could be problematic since an increasing number of studies have shown that population structure can generate spurious signatures of population size change. Here, we simulate microsatellite data from a species subdivided into social groups where reproduction occurs according to different mating systems (monogamy, polygynandry, and polygyny). We estimate the effective population size (N e) and quantify the effect of social structure on estimates of changes in N e. We analyze the simulated data with two widely used methods for demographic inference. The first approach, BOTTLENECK, tests whether the samples are at mutation-drift equilibrium and thus whether a single N e can be estimated. The second approach, msvar, aims at quantifying and dating changes in N e. We find that social structure may lead to signals of departure from mutation-drift equilibrium including signals of expansion and bottlenecks. We also find that expansion signals may be observed under simple stationary Wright-Fisher models with low diversity. Since small populations tend to characterize many endangered species, we stress that methods trying to infer N e should be interpreted with care and validated with simulated data incorporating information about structure. Spurious expansion signals due to social structure can mask critical population size changes. These can obscure true bottleneck events and be particularly problematic in endangered species.

Neutral competition explains the clonal composition of neural organoids

Neural organoids model the development of the human brain and are an indispensable tool for studying neurodevelopment. Whole-organoid lineage tracing has revealed the number of progenies arising from each initial stem cell to be highly diverse, with lineage sizes ranging from one to more than 20,000 cells. This high variability exceeds what can be explained by existing stochastic models of corticogenesis and indicates the existence of an additional source of stochasticity. To explain this variability, we introduce the SAN model which distinguishes Symmetrically diving, Asymmetrically dividing, and Non-proliferating cells. In the SAN model, the additional source of stochasticity is the survival time of a lineage's pool of symmetrically dividing cells. These survival times result from neutral competition within the sub-population of all symmetrically dividing cells. We demonstrate that our model explains the experimentally observed variability of lineage sizes and derive the quantitative relationship between survival time and lineage size. We also show that our model implies the existence of a regulatory mechanism which keeps the size of the symmetrically dividing cell population constant. Our results provide quantitative insight into the clonal composition of neural organoids and how it arises. This is relevant for many applications of neural organoids, and similar processes may occur in other developing tissues both in vitro and in vivo.

Pseudo-marginal approximation to the free energy in a micro-macro Markov chain Monte Carlo method

We introduce a generalized micro-macro Markov chain Monte Carlo (mM-MCMC) method with pseudo-marginal approximation to the free energy that is able to accelerate sampling of the microscopic Gibbs distributions when there is a time-scale separation between the macroscopic dynamics of a reaction coordinate and the remaining microscopic degrees of freedom. The mM-MCMC method attains this efficiency by iterating four steps: (i) propose a new value of the reaction coordinate, (ii) accept or reject the macroscopic sample, (iii) run a biased simulation that creates a microscopic molecular instance that lies close to the newly sampled macroscopic reaction coordinate value, and (iv) microscopic accept/reject step for the new microscopic sample. In the present paper, we eliminate the main computational bottleneck of earlier versions of this method: the necessity to have an accurate approximation of free energy. We show that the introduction of a pseudo-marginal approximation significantly reduces the computational cost of the microscopic accept/reject step while still providing unbiased samples. We illustrate the method's behavior on several molecular systems with low-dimensional reaction coordinates.

Adaptive MCMC for Bayesian Variable Selection in Generalised Linear Models and Survival Models

Developing an efficient computational scheme for high-dimensional Bayesian variable selection in generalised linear models and survival models has always been a challenging problem due to the absence of closed-form solutions to the marginal likelihood. The Reversible Jump Markov Chain Monte Carlo (RJMCMC) approach can be employed to jointly sample models and coefficients, but the effective design of the trans-dimensional jumps of RJMCMC can be challenging, making it hard to implement. Alternatively, the marginal likelihood can be derived conditional on latent variables using a data-augmentation scheme (e.g., Pólya-gamma data augmentation for logistic regression) or using other estimation methods. However, suitable data-augmentation schemes are not available for every generalised linear model and survival model, and estimating the marginal likelihood using a Laplace approximation or a correlated pseudo-marginal method can be computationally expensive. In this paper, three main contributions are presented. Firstly, we present an extended Point-wise implementation of Adaptive Random Neighbourhood Informed proposal (PARNI) to efficiently sample models directly from the marginal posterior distributions of generalised linear models and survival models. Secondly, in light of the recently proposed approximate Laplace approximation, we describe an efficient and accurate estimation method for marginal likelihood that involves adaptive parameters. Additionally, we describe a new method to adapt the algorithmic tuning parameters of the PARNI proposal by replacing Rao-Blackwellised estimates with the combination of a warm-start estimate and the ergodic average. We present numerous numerical results from simulated data and eight high-dimensional genetic mapping data-sets to showcase the efficiency of the novel PARNI proposal compared with the baseline add-delete-swap proposal.

Quantifying reliability and data deficiency in global vertebrate population trends using the Living Planet Index

Global biodiversity is facing a crisis, which must be solved through effective policies and on-the-ground conservation. But governments, NGOs, and scientists need reliable indicators to guide research, conservation actions, and policy decisions. Developing reliable indicators is challenging because the data underlying those tools is incomplete and biased. For example, the Living Planet Index tracks the changing status of global vertebrate biodiversity, but taxonomic, geographic and temporal gaps and biases are present in the aggregated data used to calculate trends. However, without a basis for real-world comparison, there is no way to directly assess an indicator's accuracy or reliability. Instead, a modelling approach can be used. We developed a model of trend reliability, using simulated datasets as stand-ins for the "real world", degraded samples as stand-ins for indicator datasets (e.g., the Living Planet Database), and a distance measure to quantify reliability by comparing partially sampled to fully sampled trends. The model revealed that the proportion of species represented in the database is not always indicative of trend reliability. Important factors are the number and length of time series, as well as their mean growth rates and variance in their growth rates, both within and between time series. We found that many trends in the Living Planet Index need more data to be considered reliable, particularly trends across the global south. In general, bird trends are the most reliable, while reptile and amphibian trends are most in need of additional data. We simulated three different solutions for reducing data deficiency, and found that collating existing data (where available) is the most efficient way to improve trend reliability, whereas revisiting previously studied populations is a quick and efficient way to improve trend reliability until new long-term studies can be completed and made available.

Solving Geophysical Inversion Problems with Intractable Likelihoods: Linearized Gaussian Approximations Versus the Correlated Pseudo-marginal Method

A geophysical Bayesian inversion problem may target the posterior distribution of geological or hydrogeological parameters given geophysical data. To account for the scatter in the petrophysical relationship linking the target parameters to the geophysical properties, this study treats the intermediate geophysical properties as latent (unobservable) variables. To perform inversion in such a latent variable model, the intractable likelihood function of the (hydro)geological parameters given the geophysical data needs to be estimated. This can be achieved by approximation with a Gaussian probability density function based on local linearization of the geophysical forward operator, thereby, accounting for the noise in the petrophysical relationship by a corresponding addition to the data covariance matrix. The new approximate method is compared against the general correlated pseudo-marginal method, which estimates the likelihood by Monte Carlo averaging over samples of the latent variable. First, the performances of the two methods are tested on a synthetic test example, in which a multivariate Gaussian porosity field is inferred using crosshole ground-penetrating radar first-arrival travel times. For this example with rather small petrophysical uncertainty, the two methods provide near-identical estimates, while an inversion that ignores petrophysical uncertainty leads to biased estimates. The results of a sensitivity analysis are then used to suggest that the linearized Gaussian approach, while attractive due to its relative computational speed, suffers from a decreasing accuracy with increasing scatter in the petrophysical relationship. The computationally more expensive correlated pseudo-marginal method performs very well even for settings with high petrophysical uncertainty.

MLNe: Simulating and Estimating Effective Size and Migration Rate from Temporal Changes in Allele Frequencies

In studies of molecular ecology, conservation biology and evolutionary biology, the current or recent effective size (Ne) of a population is frequently estimated from the marker genotype data of two or more temporally spaced samples of individuals taken from the population. Despite the developments of numerous Bayesian, likelihood and moment estimators, only a couple of them can use both temporally and spatially spaced samples of individuals to estimate jointly the effective size (Ne) of and the migration rate (m) into a population. In this note I describe new implementations of these joint estimators of Ne and m in software MLNe which runs on multiple platforms (Windows, Mac, Linux) with or without a graphical user interface (GUI), has an integrated simulation module to simulate genotype data for investigating the impacts of various factors (such as sample size and sampling interval) on estimation precision and accuracy, exploits both Message Passing Interface (MPI) and openMP for parallel computations using multiple cores and nodes to speed up analysis. The program does not require data pre-processing and accepts multiple formats of a file of original genotype data and a file of parameters as input. The GUI facilitates data and parameter inputs and produces publication-quality output graphs, while the non-GUI version of software is convenient for batch analysis of multiple datasets as in simulations. MLNe will help advance the analysis of temporal genetic marker data for estimating Ne of and m between populations, which are important parameters that will help biologists for the conservation management of natural and managed populations. MLNe can be downloaded free from the website http://www.zsl.org/science/research/software/.

Complex model calibration through emulation, a worked example for a stochastic epidemic model

Uncertainty quantification is a formal paradigm of statistical estimation that aims to account for all uncertainties inherent in the modelling process of real-world complex systems. The methods are directly applicable to stochastic models in epidemiology, however they have thus far not been widely used in this context. In this paper, we provide a tutorial on uncertainty quantification of stochastic epidemic models, aiming to facilitate the use of the uncertainty quantification paradigm for practitioners with other complex stochastic simulators of applied systems. We provide a formal workflow including the important decisions and considerations that need to be taken, and illustrate the methods over a simple stochastic epidemic model of UK SARS-CoV-2 transmission and patient outcome. We also present new approaches to visualisation of outputs from sensitivity analyses and uncertainty quantification more generally in high input and/or output dimensions.

A brief history and popularity of methods and tools used to estimate micro-evolutionary forces

Population genetics is a field of research that predates the current generations of sequencing technology. Those approaches, that were established before massively parallel sequencing methods, have been adapted to these new marker systems (in some cases involving the development of new methods) that allow genome-wide estimates of the four major micro-evolutionary forces-mutation, gene flow, genetic drift, and selection. Nevertheless, classic population genetic markers are still commonly used and a plethora of analysis methods and programs is available for these and high-throughput sequencing (HTS) data. These methods employ various and diverse theoretical and statistical frameworks, to varying degrees of success, to estimate similar evolutionary parameters making it difficult to get a concise overview across the available approaches. Presently, reviews on this topic generally focus on a particular class of methods to estimate one or two evolutionary parameters. Here, we provide a brief history of methods and a comprehensive list of available programs for estimating micro-evolutionary forces. We furthermore analyzed their usage within the research community based on popularity (citation bias) and discuss the implications of this bias for the software community. We found that a few programs received the majority of citations, with program success being independent of both the parameters estimated and the computing platform. The only deviation from a model of exponential growth in the number of citations was found for the presence of a graphical user interface (GUI). Interestingly, no relationship was found for the impact factor of the journals, when the tools were published, suggesting accessibility might be more important than visibility.

A Gibbs sampler for a class of random convex polytopes

We present a Gibbs sampler for the Dempster-Shafer (DS) approach to statistical inference for Categorical distributions. The DS framework extends the Bayesian approach, allows in particular the use of partial prior information, and yields three-valued uncertainty assessments representing probabilities "for", "against", and "don't know" about formal assertions of interest. The proposed algorithm targets the distribution of a class of random convex polytopes which encapsulate the DS inference. The sampler relies on an equivalence between the iterative constraints of the vertex configuration and the non-negativity of cycles in a fully connected directed graph. Illustrations include the testing of independence in 2 × 2 contingency tables and parameter estimation of the linkage model.

Inferring the Allelic Series at QTL in Multiparental Populations

Multiparental populations (MPPs) are experimental populations in which the genome of every individual is a mosaic of known founder haplotypes. These populations are useful for detecting quantitative trait loci (QTL) because tests of association can leverage inferred founder haplotype descent. It is difficult, however, to determine how haplotypes at a locus group into distinct functional alleles, termed the allelic series. The allelic series is important because it provides information about the number of causal variants at a QTL and their combined effects. In this study, we introduce a fully Bayesian model selection framework for inferring the allelic series. This framework accounts for sources of uncertainty found in typical MPPs, including the number and composition of functional alleles. Our prior distribution for the allelic series is based on the Chinese restaurant process, a relative of the Dirichlet process, and we leverage its connection to the coalescent to introduce additional prior information about haplotype relatedness via a phylogenetic tree. We evaluate our approach via simulation and apply it to QTL from two MPPs: the Collaborative Cross (CC) and the Drosophila Synthetic Population Resource (DSPR). We find that, although posterior inference of the exact allelic series is often uncertain, we are able to distinguish biallelic QTL from more complex multiallelic cases. Additionally, our allele-based approach improves haplotype effect estimation when the true number of functional alleles is small. Our method, Tree-Based Inference of Multiallelism via Bayesian Regression (TIMBR), provides new insight into the genetic architecture of QTL in MPPs.

Detecting and Quantifying Natural Selection at Two Linked Loci from Time Series Data of Allele Frequencies with Forward-in-Time Simulations

Recent advances in DNA sequencing techniques have made it possible to monitor genomes in great detail over time. This improvement provides an opportunity for us to study natural selection based on time serial samples of genomes while accounting for genetic recombination effect and local linkage information. Such time series genomic data allow for more accurate estimation of population genetic parameters and hypothesis testing on the recent action of natural selection. In this work, we develop a novel Bayesian statistical framework for inferring natural selection at a pair of linked loci by capitalising on the temporal aspect of DNA data with the additional flexibility of modeling the sampled chromosomes that contain unknown alleles. Our approach is built on a hidden Markov model where the underlying process is a two-locus Wright-Fisher diffusion with selection, which enables us to explicitly model genetic recombination and local linkage. The posterior probability distribution for selection coefficients is computed by applying the particle marginal Metropolis-Hastings algorithm, which allows us to efficiently calculate the likelihood. We evaluate the performance of our Bayesian inference procedure through extensive simulations, showing that our approach can deliver accurate estimates of selection coefficients, and the addition of genetic recombination and local linkage brings about significant improvement in the inference of natural selection. We also illustrate the utility of our method on real data with an application to ancient DNA data associated with white spotting patterns in horses.

Insights from Population Genomics to Enhance and Sustain Biological Control of Insect Pests

Biological control-the use of organisms (e.g., nematodes, arthropods, bacteria, fungi, viruses) for the suppression of insect pest species-is a well-established, ecologically sound and economically profitable tactic for crop protection. This approach has served as a sustainable solution for many insect pest problems for over a century in North America. However, all pest management tactics have associated risks. Specifically, the ecological non-target effects of biological control have been examined in numerous systems. In contrast, the need to understand the short- and long-term evolutionary consequences of human-mediated manipulation of biological control organisms for importation, augmentation and conservation biological control has only recently been acknowledged. Particularly, population genomics presents exceptional opportunities to study adaptive evolution and invasiveness of pests and biological control organisms. Population genomics also provides insights into (1) long-term biological consequences of releases, (2) the ecological success and sustainability of this pest management tactic and (3) non-target effects on native species, populations and ecosystems. Recent advances in genomic sequencing technology and model-based statistical methods to analyze population-scale genomic data provide a much needed impetus for biological control programs to benefit by incorporating a consideration of evolutionary consequences. Here, we review current technology and methods in population genomics and their applications to biological control and include basic guidelines for biological control researchers for implementing genomic technology and statistical modeling.

Bayesian model discrimination for partially-observed epidemic models

An efficient method for Bayesian model selection is presented for a broad class of continuous-time Markov chain models and is subsequently applied to two important problems in epidemiology. The first problem is to identify the shape of the infectious period distribution; the second problem is to determine whether individuals display symptoms before, at the same time, or after they become infectious. In both cases we show that the correct model can be identified, in the majority of cases, from symptom onset data generated from multiple outbreaks in small populations. The method works by evaluating the likelihood using a particle filter that incorporates a novel importance sampling algorithm designed for partially-observed continuous-time Markov chains. This is combined with another importance sampling method to unbiasedly estimate the model evidence. These come with estimates of precision, which allow for stopping criterion to be employed. Our method is general and can be applied to a wide range of model selection problems in biological and epidemiological systems with intractable likelihood functions.

Bayesian inference on stochastic gene transcription from flow cytometry data

Motivation

Transcription in single cells is an inherently stochastic process as mRNA levels vary greatly between cells, even for genetically identical cells under the same experimental and environmental conditions. We present a stochastic two-state switch model for the population of mRNA molecules in single cells where genes stochastically alternate between a more active ON state and a less active OFF state. We prove that the stationary solution of such a model can be written as a mixture of a Poisson and a Poisson-beta probability distribution. This finding facilitates inference for single cell expression data, observed at a single time point, from flow cytometry experiments such as FACS or fluorescence in situ hybridization (FISH) as it allows one to sample directly from the equilibrium distribution of the mRNA population. We hence propose a Bayesian inferential methodology using a pseudo-marginal approach and a recent approximation to integrate over unobserved states associated with measurement error.

Results

We provide a general inferential framework which can be widely used to study transcription in single cells from the kind of data arising in flow cytometry experiments. The approach allows us to separate between the intrinsic stochasticity of the molecular dynamics and the measurement noise. The methodology is tested in simulation studies and results are obtained for experimental multiple single cell expression data from FISH flow cytometry experiments.

Availability and implementation

All analyses were implemented in R. Source code and the experimental data are available at https://github.com/SimoneTiberi/Bayesian-inference-on-stochastic-gene-transcription-from-flow-cytometry-data.

Supplementary information

Supplementary data are available at Bioinformatics online.

Inbreeding reduces long-term growth of Alpine ibex populations

Many studies document negative inbreeding effects on individuals, and conservation efforts to preserve rare species routinely employ strategies to reduce inbreeding. Despite this, there are few clear examples in nature of inbreeding decreasing the growth rates of populations, and the extent of population-level effects of inbreeding in the wild remains controversial. Here, we take advantage of a long-term dataset of 26 reintroduced Alpine ibex (Capra ibex ibex) populations spanning nearly 100 years to show that inbreeding substantially reduced per capita population growth rates, particularly for populations in harsher environments. Populations with high average inbreeding (F ≈ 0.2) had population growth rates reduced by 71% compared with populations with no inbreeding. Our results show that inbreeding can have long-term demographic consequences even when environmental variation is large and deleterious alleles may have been purged during bottlenecks. Thus, efforts to guard against inbreeding effects in populations of endangered species have not been misplaced.

Unbiased approximations of products of expectations

We consider the problem of approximating the product of $n$ expectations with respect to a common probability distribution $\mu$. Such products routinely arise in statistics as values of the likelihood in latent variable models. Motivated by pseudo-marginal Markov chain Monte Carlo schemes, we focus on unbiased estimators of such products. The standard approach is to sample $N$ particles from $\mu$ and assign each particle to one of the expectations; this is wasteful and typically requires the number of particles to grow quadratically with the number of expectations. We propose an alternative estimator that approximates each expectation using most of the particles while preserving unbiasedness, which is computationally more efficient when the cost of simulations greatly exceeds the cost of likelihood evaluations. We carefully study the properties of our proposed estimator, showing that in latent variable contexts it needs only ${O} (n)$ particles to match the performance of the standard approach with ${O}(n^{2})$ particles. We demonstrate the procedure on two latent variable examples from approximate Bayesian computation and single-cell gene expression analysis, observing computational gains by factors of about 25 and 450, respectively.

Considerations for monitoring population trends of colonial waterbirds using the effective number of breeders and census estimates

Detecting trends in population size fluctuations is a major focus in ecology, evolution, and conservation biology. Populations of colonial waterbirds have been monitored using demographic approaches to determine annual census size (Na). We propose the addition of genetic estimates of the effective number of breeders (Nb) as indirect measures of the risk of loss of genetic diversity to improve the evaluation of demographics and increase the accuracy of trend estimates in breeding colonies. Here, we investigated which methods of the estimation of Nb are more precise under conditions of moderate genetic diversity, limited sample sizes and few microsatellite loci, as often occurs with natural populations. We used the wood stork as a model species and we offered a workflow that researchers can follow for monitoring bird breeding colonies. Our approach started with simulations using five estimators of Nb and the theoretical results were validated with empirical data collected from breeding colonies settled in the Brazilian Pantanal wetland. In parallel, we estimated census size using a corrected method based on counting active nests. Both in simulations and in natural populations, the approximate Bayesian computation (ABC) and sibship assignment (SA) methods yielded more precise estimates than the linkage disequilibrium, heterozygosity excess, and molecular coancestry methods. In particular, the ABC method performed best with few loci and small sample sizes, while the other estimators required larger sample sizes and at least 13 loci to not underestimate Nb. Moreover, according to our Nb/Na estimates (values were often ≤0.1), the wood stork colonies evaluated could be facing the loss of genetic diversity. We demonstrate that the combination of genetic and census estimates is a useful approach for monitoring natural breeding bird populations. This methodology has been recommended for populations of rare species or with a known history of population decline to support conservation efforts.

Inferring sex-specific demographic history from SNP data

The relative female and male contributions to demography are of great importance to better understand the history and dynamics of populations. While earlier studies relied on uniparental markers to investigate sex-specific questions, the increasing amount of sequence data now enables us to take advantage of tens to hundreds of thousands of independent loci from autosomes and the X chromosome. Here, we develop a novel method to estimate effective sex ratios or ESR (defined as the female proportion of the effective population) from allele count data for each branch of a rooted tree topology that summarizes the history of the populations of interest. Our method relies on Kimura’s time-dependent diffusion approximation for genetic drift, and is based on a hierarchical Bayesian model to integrate over the allele frequencies along the branches. We show via simulations that parameters are inferred robustly, even under scenarios that violate some of the model assumptions. Analyzing bovine SNP data, we infer a strongly female-biased ESR in both dairy and beef cattle, as expected from the underlying breeding scheme. Conversely, we observe a strongly male-biased ESR in early domestication times, consistent with an easier taming and management of cows, and/or introgression from wild auroch males, that would both cause a relative increase in male effective population size. In humans, analyzing a subsample of non-African populations, we find a male-biased ESR in Oceanians that may reflect complex marriage patterns in Aboriginal Australians. Because our approach relies on allele count data, it may be applied on a wide range of species.

The history of populations and their social organization is often intricate due to breeding structures, migration patterns or population bottlenecks. Estimation of the female proportion of the effective population (sex ratio) is therefore important to better understand this underlying social structure and dynamics. This question has been mainly investigated so far by comparing genetic variation of mitochondrial DNA and the Y chromosome, two uniparentally inherited markers that reflect the demographic history of females and males, respectively. To overcome the intrinsic limitations of these genetic markers, and to take advantage of the increasing amount of sequence data, we propose a new approach that uses large numbers of independent polymorphisms from autosomes and the X chromosome to estimate sex ratios, throughout the history of populations. This method allows us to confirm a strongly female-biased sex ratio in modern dairy and beef cattle breeds. Yet, we find a strongly male-biased sex ratio during domestication times, consistent with an easier taming and management of cows, and/or introgression from wild auroch males. Analyzing human data from a sample of non-African populations, we find a male bias in Oceanians, possibly indicating complex marriage patterns among Aboriginal Australian groups.