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.

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.