Population genetics in R

Allele frequencies, F-statistics, effective size and the coalescent, each built from scratch in base R and checked against a simulation with a known answer.

Population genetic data looks deceptively simple: a table of genotypes, a handful of summary statistics, a number to put in a report. The statistics are where the difficulty hides. Fst has a ceiling that moves with marker type. Ne is not a headcount and refers to a time window that depends on the estimator. A significant Tajima’s D is produced by selection and by demography alike. None of that is visible in the output of a software package.

These tutorials build each statistic from its definition in base R, then run it against simulated data whose answer is known in advance. That is the only way to see what an estimator does when its assumptions fail, and failing assumptions are the normal condition of real projects.

Allele frequencies and population structure

The starting point is one locus in one population, and the first lesson is that the standard test is a smoke alarm rather than a diagnosis.

• Hardy-Weinberg expectations builds the genotype expectation from gene counting, implements the exact test that replaces the chi-square when an allele is rare, and shows two unrelated mechanisms producing the same heterozygote deficit.
• F-statistics and population structure derives Fis, Fst and Fit from three heterozygosities, demonstrates why a ratio of heterozygosities cannot compare marker types, and swaps the descriptive ratio for an estimator that survives small samples.
• Drift, migration and isolation by distance simulates both forces, confirms the classical migration-drift equilibrium, and works out how little a measured Fst licenses you to say about the number of migrants.
• Checking a population genetics analysis is the diagnostic set: how much information the design carries, how to tell a marker artefact from real structure, what a routine filtering step does to the estimate, and when an outlier locus means nothing.

Effective size and bottlenecks

The rate at which a population loses diversity is set by its breeding arrangement rather than its census, and estimating that rate from markers is a separate problem with its own failure modes.

• Effective population size measures the loss of diversity in simulated populations of identical census size, then checks the textbook formulas for sex ratio, reproductive variance and fluctuating size one at a time.
• Estimating Ne from genetic data implements the temporal and linkage-disequilibrium estimators, calibrates the sampling floor by permutation, and finds where each method stops working.
• Bottlenecks and genetic diversity pushes a population through crashes of different depths and compares what allelic richness, heterozygosity and the allele size range notice.
• Checking an effective size estimate asks which time window an estimate covers, what sampling siblings does to it, how wide the interval really is, and what physical linkage costs.

Coalescent theory and sequence diversity

Run the same process backwards and the ancestry of the sample becomes the object of study. Sequence data is mutations on that unobserved tree.

• The coalescent and gene trees builds the backward process from exponential waiting times, verifies tree height and total branch length, and shows that a single gene tree is one draw from a very broad distribution.
• Nucleotide diversity and Watterson theta places mutations on the branches, derives the two classical estimators of the scaled mutation parameter, and explains which part of the tree each one trusts.
• Tajima’s D and the frequency spectrum turns the disagreement between those estimators into a statistic, then breaks the neutral model with growth, a bottleneck and population structure to see which direction each one pushes.
• Checking a coalescent analysis measures what one locus can support, what recombination buys for free, whether demography really imitates selection, and where to spend the sequencing budget.

Clustering, ancestry and assignment

The bar plot that ends up in almost every population genetics paper is a model fit, and it is worth knowing which parts of it are estimates and which are artefacts of the algorithm.

• Admixture and ancestry proportions in R fits the K-source ancestry model by expectation maximisation and measures the three things that decide whether the resulting bar plot carries information: the number of loci, how differentiated the sources are, and the bias that keeps a genuinely unmixed individual off the boundary.
• Choosing K in genetic clustering compares the likelihood, cross-validation and the delta-K heuristic on data with real clusters and then on a cline that has none, where all three still recommend more than one group.
• Assignment tests and self-assignment asks the narrower question that has a measurable error rate, separates honest leave-one-out accuracy from the flattering version, and handles the source population nobody sampled.
• Checking a genetic clustering analysis works through uneven sample sizes, linked loci that count twice, the spread across random starts, and what survives dropping half the loci or a whole source.

What these tutorials assume

Base R and ggplot2, nothing else. Every simulation is written out rather than called from a package, because the point is to see the assumption that the package hides. The data is simulated throughout, so every estimate can be compared with the truth, and the numbers quoted in the prose are the numbers the code produces when you run it.

Two conventions are worth stating. Diploid, biallelic markers are the default, with multi-allelic and stepwise-mutation models introduced where they matter. And time is measured in generations, or in units of twice the effective size where the coalescent makes that natural, which is the same statement in different currency.

Where this connects

The methods here overlap several other clusters on this site, and the overlaps are worth following rather than avoiding. The distance-matrix machinery behind isolation-by-distance analyses lives with the Mantel test tutorials in the diversity cluster, which is where the caveats about correlating two distance matrices belong.

• Quantitative genetics in R works on trait variation rather than allele frequencies: selection differentials and gradients, the breeder’s equation, and the genetic constraints that decide which direction a population can respond in.
• Relatedness as a regression treats relatedness as a statistical coefficient rather than a pedigree fraction, which is the same move that turns Fis into a correlation between uniting gametes.
• The Moran process and fixation probability and its cluster follow the same binomial sampling forwards in small populations, where the outcome of interest is fixation rather than differentiation.
• Landscape connectivity in R sends the same differentiation estimates outdoors, regressing them on straight-line, least-cost and resistance distances and asking which map of movement cost the genetic data support.