Checking a population genetics analysis

population genetics

model diagnostics

ecology tutorial

ggplot2

Four checks for an Fst analysis in R: loci against individuals, separating substructure from marker error, filtering rare variants, and outlier tests.

Author

Tidy Ecology

Published

2026-08-11

The three tutorials before this one worked with data whose answer was known. That is the only way to learn what an estimator does, and it is also why the numbers came out reassuring. This one is about the checks that stay useful when the answer is not known, which is every real project.

Four checks, in the order they usually matter: whether the design has enough information, whether a positive Fis is biology or a marker problem, what a routine filtering step does to the estimate, and whether an outlier locus means selection. Each one is a simulation with a known truth, run so that the diagnostic can be calibrated before it is trusted.

library(ggplot2)

te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
               clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
               ink = "#16241d")

theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e6dc"),
          plot.background = element_rect(fill = "white", colour = NA),
          panel.background = element_rect(fill = "white", colour = NA),
          plot.title = element_text(face = "bold", colour = te_pal$ink),
          axis.title = element_text(colour = "#2c3a31"))
}

Two helpers carry the whole post. The first simulates many biallelic loci in several populations with a known Fst, using the beta model from the second tutorial and drawing genotype counts directly. The second is the Weir and Cockerham estimator, vectorised across loci, combining them as a ratio of averages.

sim_loci <- function(n_loci, r, n_ind, fst, p_anc) {
  if (length(p_anc) == 1) p_anc <- rep(p_anc, n_loci)
  aa <- ab <- bb <- matrix(0L, n_loci, r)
  shape1 <- p_anc * (1 - fst) / fst
  shape2 <- (1 - p_anc) * (1 - fst) / fst
  for (j in seq_len(r)) {
    p_j <- rbeta(n_loci, shape1, shape2)
    het <- rbinom(n_loci, n_ind, 2 * p_j * (1 - p_j))
    rest <- n_ind - het
    hom_a <- rbinom(n_loci, rest, p_j^2 / (p_j^2 + (1 - p_j)^2))
    aa[, j] <- hom_a
    ab[, j] <- het
    bb[, j] <- rest - hom_a
  }
  list(aa = aa, ab = ab, bb = bb, n = n_ind, r = r)
}

wc_multi <- function(g) {
  n <- g$n
  r <- g$r
  p <- (2 * g$aa + g$ab) / (2 * n)
  h <- g$ab / n
  p_w <- rowMeans(p)
  s2 <- rowSums((p - p_w)^2) / (r - 1)
  h_w <- rowMeans(h)
  a <- s2 - (p_w * (1 - p_w) - (r - 1) / r * s2 - h_w / 4) / (n - 1)
  b <- n / (n - 1) * (p_w * (1 - p_w) - (r - 1) / r * s2 -
                        (2 * n - 1) / (4 * n) * h_w)
  cw <- h_w / 2
  list(theta = sum(a) / sum(a + b + cw), per_locus = a / (a + b + cw))
}

Check 1: loci or individuals?

Field time buys individuals and laboratory budget buys loci, and the two are not interchangeable. Fix the truth at Fst = 0.05 in five populations and measure the precision of the multilocus estimate across a grid.

set.seed(471)

n_rep <- 200
design <- expand.grid(n_ind = c(2, 4, 8, 16, 32), n_loci = c(10, 100, 1000))
design$mean_theta <- NA_real_
design$sd_theta <- NA_real_
for (k in seq_len(nrow(design))) {
  est <- replicate(n_rep, wc_multi(sim_loci(design$n_loci[k], 5,
                                            design$n_ind[k], 0.05, 0.3))$theta)
  design$mean_theta[k] <- mean(est)
  design$sd_theta[k] <- sd(est)
}
design$sd_times_root_nL <- design$sd_theta * sqrt(design$n_ind * design$n_loci)
print(round(design, 4))

   n_ind n_loci mean_theta sd_theta sd_times_root_nL
1      2     10     0.0547   0.0763           0.3413
2      4     10     0.0464   0.0396           0.2505
3      8     10     0.0499   0.0242           0.2162
4     16     10     0.0493   0.0170           0.2156
5     32     10     0.0502   0.0130           0.2333
6      2    100     0.0509   0.0257           0.3634
7      4    100     0.0517   0.0133           0.2657
8      8    100     0.0504   0.0078           0.2212
9     16    100     0.0501   0.0053           0.2133
10    32    100     0.0499   0.0043           0.2405
11     2   1000     0.0489   0.0071           0.3164
12     4   1000     0.0497   0.0035           0.2245
13     8   1000     0.0498   0.0024           0.2185
14    16   1000     0.0500   0.0016           0.2011
15    32   1000     0.0501   0.0013           0.2395

Two things to read off. First, the mean estimate is centred on the truth in every cell, from 0.0464 to 0.0547 across the whole grid and from 0.0489 to 0.0517 once there are 100 loci, including the cells with two individuals per population: the estimator is unbiased regardless of sample size, which is the property Nei’s ratio lacked. Only the scatter changes.

Second, the scatter is governed by the product. The last column is the standard deviation multiplied by the square root of the number of genotypes, and from four individuals upward it sits between 0.2011 and 0.2657 everywhere in the grid, across designs whose smallest and largest differ by nearly three orders of magnitude in size. Precision therefore comes from n * L, and the design question becomes which factor is cheaper to raise. Two individuals do pay a penalty (0.3164 to 0.3634), so a floor of about four individuals per population is worth respecting, but above that floor it makes no difference to the arithmetic whether the genotypes come from many individuals or many loci.

ggplot(design, aes(n_loci, sd_theta, colour = factor(n_ind))) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 2) +
  scale_x_log10(breaks = c(10, 100, 1000)) +
  scale_y_log10() +
  scale_colour_manual(values = c(`2` = te_pal$clay, `4` = te_pal$gold,
                                 `8` = te_pal$sage, `16` = te_pal$green,
                                 `32` = te_pal$forest),
                      name = "individuals\nper population") +
  labs(title = "Precision follows the number of genotypes",
       subtitle = "true Fst 0.05, five populations, 200 replicates per point",
       x = "loci (log scale)", y = "SD of the estimate (log scale)") +
  theme_te()

Five parallel downward lines on log axes, one per sample size, showing the standard deviation of Fst falling from about 0.08 to 0.001 as loci increase from ten to a thousand. — Figure 1: Standard deviation of the multilocus Fst estimate against the number of loci, for five sample sizes, on log axes. True Fst is 0.05 in five populations, 200 replicates per point.

The caveat is the part the simulation cannot show. Here every locus receives its own independent draw of population history, so adding loci keeps shrinking the error without limit. Real loci share one history and are linked into blocks, so beyond a few thousand markers the error stops falling: what remains is the uncertainty of the single realised history, and no amount of sequencing removes it. Treat the grid as an upper bound on what loci buy.

Check 2: is that Fis biology or a pipette?

A positive Fis was ambiguous in the first tutorial: population structure and allelic dropout produce the same heterozygote deficit. With more than one locus the ambiguity is resolvable, because the two mechanisms leave different footprints across loci. Structure shifts every locus a little; a marker artefact wrecks a few loci badly.

Both scenarios below are built to give almost the same mean Fis, which is what makes the point.

fis_per_locus <- function(aa, ab, bb) {
  n <- aa + ab + bb
  p <- (2 * aa + ab) / (2 * n)
  1 - (ab / n) / (2 * p * (1 - p))
}

set.seed(472)
n_loci <- 100
n_ind <- 400

## scenario A: one sample that is really two subpopulations, Fst 0.015
fst_sub <- 0.015
p_anc <- runif(n_loci, 0.15, 0.5)
sh1 <- p_anc * (1 - fst_sub) / fst_sub
sh2 <- (1 - p_anc) * (1 - fst_sub) / fst_sub
p1 <- rbeta(n_loci, sh1, sh2)
p2 <- rbeta(n_loci, sh1, sh2)

mix_geno <- function(p_a, p_b, n_half) {
  het <- rbinom(length(p_a), n_half, 2 * p_a * (1 - p_a)) +
    rbinom(length(p_b), n_half, 2 * p_b * (1 - p_b))
  hom_a <- rbinom(length(p_a), n_half, p_a^2) +
    rbinom(length(p_b), n_half, p_b^2)
  cbind(aa = hom_a, ab = het, bb = 2 * n_half - hom_a - het)
}

sub_dat <- mix_geno(p1, p2, n_ind / 2)
fis_sub <- fis_per_locus(sub_dat[, "aa"], sub_dat[, "ab"], sub_dat[, "bb"])

## scenario B: random mating, but five loci lose 30 percent of heterozygotes
het <- rbinom(n_loci, n_ind, 2 * p_anc * (1 - p_anc))
hom_a <- rbinom(n_loci, n_ind - het, p_anc^2 / (p_anc^2 + (1 - p_anc)^2))
hom_b <- n_ind - het - hom_a
bad <- 1:5
lost <- rbinom(length(bad), het[bad], 0.3)
lost_a <- rbinom(length(bad), lost, 0.5)
het[bad] <- het[bad] - lost
hom_a[bad] <- hom_a[bad] + lost_a
hom_b[bad] <- hom_b[bad] + lost - lost_a
fis_art <- fis_per_locus(hom_a, het, hom_b)

report <- function(x, label) {
  cat(label, " mean", round(mean(x), 4), " median", round(median(x), 4),
      " sd", round(sd(x), 4), " max", round(max(x), 4),
      " loci above 0.15:", sum(x > 0.15),
      " loci above zero:", sum(x > 0), "\n")
}
report(fis_sub, "hidden substructure:")

hidden substructure:  mean 0.0107  median 0.0156  sd 0.0515  max 0.1146  loci above 0.15: 0  loci above zero: 59

report(fis_art, "dropout at 5 loci: ")

dropout at 5 loci:   mean 0.0123  median -0.005  sd 0.0881  max 0.3824  loci above 0.15: 5  loci above zero: 48

cat("the five damaged loci:", round(fis_art[bad], 3), "\n")

the five damaged loci: 0.364 0.295 0.382 0.211 0.332

cat("the other 95 loci: mean", round(mean(fis_art[-bad]), 4),
    " median", round(median(fis_art[-bad]), 4), "\n")

the other 95 loci: mean -0.0037  median -0.0089

Mean Fis is 0.0107 under hidden substructure and 0.0123 under dropout, so a report that quotes only the mean cannot distinguish them. Everything else does. Under substructure the median is 0.0156, positive like the mean, 59 of 100 loci are above zero, and the largest single value is 0.1146. Under dropout the median is -0.005, the undamaged 95 loci average -0.0037, and the five broken loci come in at 0.211 to 0.382.

fis_df <- data.frame(
  locus = rep(1:n_loci, 2),
  fis = c(fis_sub, fis_art),
  scenario = rep(c("hidden substructure", "dropout at 5 loci"),
                 each = n_loci))
fis_df$scenario <- factor(fis_df$scenario,
                          levels = c("hidden substructure",
                                     "dropout at 5 loci"))

ggplot(fis_df, aes(locus, fis)) +
  geom_hline(yintercept = 0, colour = te_pal$line) +
  geom_point(aes(colour = scenario), size = 1.7, show.legend = FALSE) +
  facet_wrap(~scenario) +
  scale_colour_manual(values = c(`hidden substructure` = te_pal$forest,
                                 `dropout at 5 loci` = te_pal$clay)) +
  labs(title = "Same mean Fis, different footprint across loci",
       subtitle = "100 loci, 400 individuals",
       x = "locus", y = "per-locus Fis") +
  theme_te()

Two panels of per-locus points around zero; the substructure panel is shifted slightly upward as a whole, the dropout panel is centred on zero with five points far above the rest. — Figure 2: Per-locus Fis under hidden substructure and under allelic dropout at five loci. The two data sets have nearly the same mean Fis.

The check is therefore to plot Fis per locus and never to average first. A shifted cloud says the sample is not one randomly mating population, and the remedy is to look for the grouping. A flat cloud with a few extreme loci says those loci are broken, and the remedy is to drop them. A biologically inbred population gives the shifted cloud too, so this check narrows the possibilities to two rather than one; separating true inbreeding from unrecognised structure needs the sampling design, not the genotypes.

Check 3: what filtering rare variants does

Almost every pipeline discards markers below a minor allele frequency threshold. That step is not neutral for Fst, because rare variants and common variants carry different amounts of among-population variance. Simulate a realistic frequency spectrum with many rare variants, then estimate Fst on the whole set and on the filtered set.

set.seed(473)

p_spec <- 1 / seq(1, 40)
p_spec <- sample(p_spec[p_spec >= 0.01], 4000, replace = TRUE)
g_all <- sim_loci(4000, 5, 25, 0.05, p_spec)

p_pooled <- rowMeans((2 * g_all$aa + g_all$ab) / (2 * g_all$n))
maf <- pmin(p_pooled, 1 - p_pooled)
keep <- maf >= 0.05

subset_g <- function(g, idx) {
  list(aa = g$aa[idx, ], ab = g$ab[idx, ], bb = g$bb[idx, ], n = g$n, r = g$r)
}

cat("all loci:", nrow(g_all$aa),
    " theta:", round(wc_multi(g_all)$theta, 4), "\n")

all loci: 4000  theta: 0.0503

cat("loci kept at MAF 0.05:", sum(keep),
    " theta:", round(wc_multi(subset_g(g_all, keep))$theta, 4), "\n")

loci kept at MAF 0.05: 1964  theta: 0.0547

cat("loci discarded:", sum(!keep),
    " theta:", round(wc_multi(subset_g(g_all, !keep))$theta, 4), "\n")

loci discarded: 2036  theta: 0.0309

The truth is 0.05 and the full set returns 0.0503. Filtering at a minor allele frequency of 0.05 throws away 2036 of 4000 loci and returns 0.0547, nine percent high, while the discarded loci on their own give 0.0309. The mechanism is not mysterious: a locus whose minor allele is rare has little variance to divide up, and the sampling noise in its frequency estimate is a larger share of what is left. Filtering removes those loci and the average shifts.

The consequence is a comparability rule. An Fst from a filtered SNP set is not on the same scale as an Fst from an unfiltered one, and a threshold applied to the pooled sample is not the same as a threshold applied per population. When two studies of the same species disagree about Fst, the filtering step is one of the first places to look, alongside the marker type argument from the second tutorial.

Check 4: an outlier is only an outlier under a model

Fst outlier scans look for loci more differentiated than a neutral null and call them candidates for selection. The scan is only as good as the null, and the usual null is an island model in which every population is equally related to every other. Real sampling rarely looks like that: populations fall into regions, catchments, or lineages.

Simulate two groups of three demes, with real structure among the groups and very little within them, and no selection anywhere. Then run the standard scan with an island-model null calibrated to the same overall Fst.

set.seed(474)

hier <- function(n_loci, n_ind, fst_among, fst_within, p_anc) {
  sh1 <- p_anc * (1 - fst_among) / fst_among
  sh2 <- (1 - p_anc) * (1 - fst_among) / fst_among
  out_aa <- out_ab <- out_bb <- matrix(0L, n_loci, 6)
  for (grp in 1:2) {
    p_grp <- rbeta(n_loci, sh1, sh2)
    s1 <- p_grp * (1 - fst_within) / fst_within
    s2 <- (1 - p_grp) * (1 - fst_within) / fst_within
    for (d in 1:3) {
      col <- (grp - 1) * 3 + d
      p_d <- rbeta(n_loci, s1, s2)
      het <- rbinom(n_loci, n_ind, 2 * p_d * (1 - p_d))
      rest <- n_ind - het
      hom_a <- rbinom(n_loci, rest, p_d^2 / (p_d^2 + (1 - p_d)^2))
      out_aa[, col] <- hom_a
      out_ab[, col] <- het
      out_bb[, col] <- rest - hom_a
    }
  }
  list(aa = out_aa, ab = out_ab, bb = out_bb, n = n_ind, r = 6)
}

p_base <- runif(2000, 0.1, 0.5)
obs <- hier(2000, 30, 0.08, 0.01, p_base)
obs_fit <- wc_multi(obs)
cat("hierarchical data, multilocus theta:", round(obs_fit$theta, 4), "\n")

hierarchical data, multilocus theta: 0.0592

null_g <- sim_loci(20000, 6, 30, obs_fit$theta, sample(p_base, 20000, TRUE))
null_per <- wc_multi(null_g)$per_locus
thr <- quantile(null_per, 0.99, na.rm = TRUE)

cat("island-model null, 99th percentile of per-locus theta:",
    round(thr, 4), "\n")

island-model null, 99th percentile of per-locus theta: 0.1949

cat("hierarchical loci above that threshold:",
    round(mean(obs_fit$per_locus > thr), 4),
    " nominal rate: 0.01\n")

hierarchical loci above that threshold: 0.048  nominal rate: 0.01

cat("null spread: sd", round(sd(null_per), 4),
    " hierarchical spread: sd", round(sd(obs_fit$per_locus), 4), "\n")

null spread: sd 0.0445  hierarchical spread: sd 0.067

Not one locus in this simulation is under selection. The scan flags 4.8 percent of them, nearly five times the nominal one percent, because the hierarchical arrangement widens the distribution of per-locus Fst beyond what the island model expects (standard deviation 0.067 against 0.0445). Every flagged locus would be reported as a candidate, and a gene ontology analysis of the list would find something to say about it.

dens_df <- data.frame(
  theta = c(null_per, obs_fit$per_locus),
  model = rep(c("island-model null", "hierarchical truth"),
              c(length(null_per), length(obs_fit$per_locus))))

ggplot(dens_df, aes(theta, colour = model)) +
  geom_density(linewidth = 0.9) +
  geom_vline(xintercept = thr, linetype = "dashed", colour = "grey40") +
  coord_cartesian(xlim = c(-0.05, 0.45)) +
  scale_colour_manual(values = c(`island-model null` = te_pal$sage,
                                 `hierarchical truth` = te_pal$clay),
                      name = NULL) +
  labs(title = "The same overall Fst, a heavier tail",
       subtitle = "dashed line: 99th percentile of the island-model null",
       x = "per-locus Weir and Cockerham estimate", y = "density") +
  theme_te()

Two overlapping density curves; the hierarchical curve has a heavier right tail than the island-model null, and a dashed vertical line at the null's 99th percentile cuts through it. — Figure 3: Per-locus Fst under an island-model null and under a hierarchical structure with the same overall Fst, with the 99th percentile of the null marked.

The fix is a null that matches the sampling: a hierarchical island model with the groups declared, or a demographic model fitted to the data first, or a scan that conditions on the observed covariance among populations. None of these is a formality. An outlier scan is a comparison against a model of the population’s history, so it answers “is this locus unusual given that history”, and if the history is wrong the answer is about the history.

The honest limit

The four checks constrain a population genetic analysis; they do not validate it. Two things they cannot reach are worth stating plainly.

The first is that no property of allele frequencies separates selection from demography at a single locus. A bottleneck, a range expansion, or an unsampled ghost population all produce loci with extreme Fst, and they do it without any locus being special. The evidence that promotes a candidate from a list to a result comes from outside the scan: a repeated signal in an independent population pair, a match to an environmental gradient, a function that makes sense, an experiment.

The second is that all four checks assume the marker set is a fair sample of the genome. Markers chosen because they were variable in one population, or an assay developed on one lineage, tilt every statistic in this cluster in the same direction, and no internal diagnostic can see it, because the data set does not contain the loci that were never assayed. That is the same shape of problem as the omitted trait in a selection analysis: the fix is in the design, not in the diagnostics.

Where to go next

That closes the population genetics cluster: expectations for one locus in one population, F-statistics for many populations, the drift and migration processes that predict Fst, and the checks. The natural continuation is the other place where allele frequencies change on their own, the stochastic dynamics of finite populations, where the same binomial sampling that produced drift here drives fixation probabilities and evolutionary outcomes.

References

Weir BS, Cockerham CC 1984 Evolution 38(6):1358-1370 (10.1111/j.1558-5646.1984.tb05657.x)

Beaumont MA, Nichols RA 1996 Proceedings of the Royal Society B 263(1377):1619-1626 (10.1098/rspb.1996.0237)

Excoffier L, Hofer T, Foll M 2009 Heredity 103(4):285-298 (10.1038/hdy.2009.74)

Willing EM, Dreyer C, van Oosterhout C 2012 PLoS ONE 7(8):e42649 (10.1371/journal.pone.0042649)

Lotterhos KE, Whitlock MC 2014 Molecular Ecology 23(9):2178-2192 (10.1111/mec.12725)

Bhatia G, Patterson N, Sankararaman S, Price AL 2013 Genome Research 23(9):1514-1521 (10.1101/gr.154831.113)

Dakin EE, Avise JC 2004 Heredity 93(5):504-509 (10.1038/sj.hdy.6800545)

--- title: "Checking a population genetics analysis" description: "Four checks for an Fst analysis in R: loci against individuals, separating substructure from marker error, filtering rare variants, and outlier tests." date: "2026-08-11 12:00" categories: [R, population genetics, model diagnostics, ecology tutorial, ggplot2] image: thumbnail.png image-alt: "Per-locus Fis values for two data sets with the same mean, one shifted as a whole and one flat with five extreme outliers." --- The three tutorials before this one worked with data whose answer was known. That is the only way to learn what an estimator does, and it is also why the numbers came out reassuring. This one is about the checks that stay useful when the answer is not known, which is every real project. Four checks, in the order they usually matter: whether the design has enough information, whether a positive `Fis` is biology or a marker problem, what a routine filtering step does to the estimate, and whether an outlier locus means selection. Each one is a simulation with a known truth, run so that the diagnostic can be calibrated before it is trusted. ```{r setup} #| message: false #| warning: false library(ggplot2) te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f", clay = "#b5534e", gold = "#cda23f", line = "#dad9ca", ink = "#16241d") theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e6dc"), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA), plot.title = element_text(face = "bold", colour = te_pal$ink), axis.title = element_text(colour = "#2c3a31")) } ``` Two helpers carry the whole post. The first simulates many biallelic loci in several populations with a known `Fst`, using the beta model from the second tutorial and drawing genotype counts directly. The second is the Weir and Cockerham estimator, vectorised across loci, combining them as a ratio of averages. ```{r helpers} sim_loci <- function(n_loci, r, n_ind, fst, p_anc) { if (length(p_anc) == 1) p_anc <- rep(p_anc, n_loci) aa <- ab <- bb <- matrix(0L, n_loci, r) shape1 <- p_anc * (1 - fst) / fst shape2 <- (1 - p_anc) * (1 - fst) / fst for (j in seq_len(r)) { p_j <- rbeta(n_loci, shape1, shape2) het <- rbinom(n_loci, n_ind, 2 * p_j * (1 - p_j)) rest <- n_ind - het hom_a <- rbinom(n_loci, rest, p_j^2 / (p_j^2 + (1 - p_j)^2)) aa[, j] <- hom_a ab[, j] <- het bb[, j] <- rest - hom_a } list(aa = aa, ab = ab, bb = bb, n = n_ind, r = r) } wc_multi <- function(g) { n <- g$n r <- g$r p <- (2 * g$aa + g$ab) / (2 * n) h <- g$ab / n p_w <- rowMeans(p) s2 <- rowSums((p - p_w)^2) / (r - 1) h_w <- rowMeans(h) a <- s2 - (p_w * (1 - p_w) - (r - 1) / r * s2 - h_w / 4) / (n - 1) b <- n / (n - 1) * (p_w * (1 - p_w) - (r - 1) / r * s2 - (2 * n - 1) / (4 * n) * h_w) cw <- h_w / 2 list(theta = sum(a) / sum(a + b + cw), per_locus = a / (a + b + cw)) } ``` ## Check 1: loci or individuals? Field time buys individuals and laboratory budget buys loci, and the two are not interchangeable. Fix the truth at `Fst` = 0.05 in five populations and measure the precision of the multilocus estimate across a grid. ```{r design} set.seed(471) n_rep <- 200 design <- expand.grid(n_ind = c(2, 4, 8, 16, 32), n_loci = c(10, 100, 1000)) design$mean_theta <- NA_real_ design$sd_theta <- NA_real_ for (k in seq_len(nrow(design))) { est <- replicate(n_rep, wc_multi(sim_loci(design$n_loci[k], 5, design$n_ind[k], 0.05, 0.3))$theta) design$mean_theta[k] <- mean(est) design$sd_theta[k] <- sd(est) } design$sd_times_root_nL <- design$sd_theta * sqrt(design$n_ind * design$n_loci) print(round(design, 4)) ``` Two things to read off. First, the mean estimate is centred on the truth in every cell, from 0.0464 to 0.0547 across the whole grid and from 0.0489 to 0.0517 once there are 100 loci, including the cells with two individuals per population: the estimator is unbiased regardless of sample size, which is the property Nei's ratio lacked. Only the scatter changes. Second, the scatter is governed by the product. The last column is the standard deviation multiplied by the square root of the number of genotypes, and from four individuals upward it sits between 0.2011 and 0.2657 everywhere in the grid, across designs whose smallest and largest differ by nearly three orders of magnitude in size. Precision therefore comes from `n * L`, and the design question becomes which factor is cheaper to raise. Two individuals do pay a penalty (0.3164 to 0.3634), so a floor of about four individuals per population is worth respecting, but above that floor it makes no difference to the arithmetic whether the genotypes come from many individuals or many loci. ```{r fig-precision} #| fig-width: 7 #| fig-height: 4.2 #| fig-cap: "Standard deviation of the multilocus Fst estimate against the number of loci, for five sample sizes, on log axes. True Fst is 0.05 in five populations, 200 replicates per point." #| fig-alt: "Five parallel downward lines on log axes, one per sample size, showing the standard deviation of Fst falling from about 0.08 to 0.001 as loci increase from ten to a thousand." ggplot(design, aes(n_loci, sd_theta, colour = factor(n_ind))) + geom_line(linewidth = 0.9) + geom_point(size = 2) + scale_x_log10(breaks = c(10, 100, 1000)) + scale_y_log10() + scale_colour_manual(values = c(`2` = te_pal$clay, `4` = te_pal$gold, `8` = te_pal$sage, `16` = te_pal$green, `32` = te_pal$forest), name = "individuals\nper population") + labs(title = "Precision follows the number of genotypes", subtitle = "true Fst 0.05, five populations, 200 replicates per point", x = "loci (log scale)", y = "SD of the estimate (log scale)") + theme_te() ``` The caveat is the part the simulation cannot show. Here every locus receives its own independent draw of population history, so adding loci keeps shrinking the error without limit. Real loci share one history and are linked into blocks, so beyond a few thousand markers the error stops falling: what remains is the uncertainty of the single realised history, and no amount of sequencing removes it. Treat the grid as an upper bound on what loci buy. ## Check 2: is that Fis biology or a pipette? A positive `Fis` was ambiguous in the first tutorial: population structure and allelic dropout produce the same heterozygote deficit. With more than one locus the ambiguity is resolvable, because the two mechanisms leave different footprints across loci. Structure shifts every locus a little; a marker artefact wrecks a few loci badly. Both scenarios below are built to give almost the same mean `Fis`, which is what makes the point. ```{r fis-diagnosis} fis_per_locus <- function(aa, ab, bb) { n <- aa + ab + bb p <- (2 * aa + ab) / (2 * n) 1 - (ab / n) / (2 * p * (1 - p)) } set.seed(472) n_loci <- 100 n_ind <- 400 ## scenario A: one sample that is really two subpopulations, Fst 0.015 fst_sub <- 0.015 p_anc <- runif(n_loci, 0.15, 0.5) sh1 <- p_anc * (1 - fst_sub) / fst_sub sh2 <- (1 - p_anc) * (1 - fst_sub) / fst_sub p1 <- rbeta(n_loci, sh1, sh2) p2 <- rbeta(n_loci, sh1, sh2) mix_geno <- function(p_a, p_b, n_half) { het <- rbinom(length(p_a), n_half, 2 * p_a * (1 - p_a)) + rbinom(length(p_b), n_half, 2 * p_b * (1 - p_b)) hom_a <- rbinom(length(p_a), n_half, p_a^2) + rbinom(length(p_b), n_half, p_b^2) cbind(aa = hom_a, ab = het, bb = 2 * n_half - hom_a - het) } sub_dat <- mix_geno(p1, p2, n_ind / 2) fis_sub <- fis_per_locus(sub_dat[, "aa"], sub_dat[, "ab"], sub_dat[, "bb"]) ## scenario B: random mating, but five loci lose 30 percent of heterozygotes het <- rbinom(n_loci, n_ind, 2 * p_anc * (1 - p_anc)) hom_a <- rbinom(n_loci, n_ind - het, p_anc^2 / (p_anc^2 + (1 - p_anc)^2)) hom_b <- n_ind - het - hom_a bad <- 1:5 lost <- rbinom(length(bad), het[bad], 0.3) lost_a <- rbinom(length(bad), lost, 0.5) het[bad] <- het[bad] - lost hom_a[bad] <- hom_a[bad] + lost_a hom_b[bad] <- hom_b[bad] + lost - lost_a fis_art <- fis_per_locus(hom_a, het, hom_b) report <- function(x, label) { cat(label, " mean", round(mean(x), 4), " median", round(median(x), 4), " sd", round(sd(x), 4), " max", round(max(x), 4), " loci above 0.15:", sum(x > 0.15), " loci above zero:", sum(x > 0), "\n") } report(fis_sub, "hidden substructure:") report(fis_art, "dropout at 5 loci: ") cat("the five damaged loci:", round(fis_art[bad], 3), "\n") cat("the other 95 loci: mean", round(mean(fis_art[-bad]), 4), " median", round(median(fis_art[-bad]), 4), "\n") ``` Mean `Fis` is 0.0107 under hidden substructure and 0.0123 under dropout, so a report that quotes only the mean cannot distinguish them. Everything else does. Under substructure the median is 0.0156, positive like the mean, 59 of 100 loci are above zero, and the largest single value is 0.1146. Under dropout the median is -0.005, the undamaged 95 loci average -0.0037, and the five broken loci come in at 0.211 to 0.382. ```{r fig-fis} #| fig-width: 7 #| fig-height: 4.2 #| fig-cap: "Per-locus Fis under hidden substructure and under allelic dropout at five loci. The two data sets have nearly the same mean Fis." #| fig-alt: "Two panels of per-locus points around zero; the substructure panel is shifted slightly upward as a whole, the dropout panel is centred on zero with five points far above the rest." fis_df <- data.frame( locus = rep(1:n_loci, 2), fis = c(fis_sub, fis_art), scenario = rep(c("hidden substructure", "dropout at 5 loci"), each = n_loci)) fis_df$scenario <- factor(fis_df$scenario, levels = c("hidden substructure", "dropout at 5 loci")) ggplot(fis_df, aes(locus, fis)) + geom_hline(yintercept = 0, colour = te_pal$line) + geom_point(aes(colour = scenario), size = 1.7, show.legend = FALSE) + facet_wrap(~scenario) + scale_colour_manual(values = c(`hidden substructure` = te_pal$forest, `dropout at 5 loci` = te_pal$clay)) + labs(title = "Same mean Fis, different footprint across loci", subtitle = "100 loci, 400 individuals", x = "locus", y = "per-locus Fis") + theme_te() ``` The check is therefore to plot `Fis` per locus and never to average first. A shifted cloud says the sample is not one randomly mating population, and the remedy is to look for the grouping. A flat cloud with a few extreme loci says those loci are broken, and the remedy is to drop them. A biologically inbred population gives the shifted cloud too, so this check narrows the possibilities to two rather than one; separating true inbreeding from unrecognised structure needs the sampling design, not the genotypes. ## Check 3: what filtering rare variants does Almost every pipeline discards markers below a minor allele frequency threshold. That step is not neutral for `Fst`, because rare variants and common variants carry different amounts of among-population variance. Simulate a realistic frequency spectrum with many rare variants, then estimate `Fst` on the whole set and on the filtered set. ```{r maf} set.seed(473) p_spec <- 1 / seq(1, 40) p_spec <- sample(p_spec[p_spec >= 0.01], 4000, replace = TRUE) g_all <- sim_loci(4000, 5, 25, 0.05, p_spec) p_pooled <- rowMeans((2 * g_all$aa + g_all$ab) / (2 * g_all$n)) maf <- pmin(p_pooled, 1 - p_pooled) keep <- maf >= 0.05 subset_g <- function(g, idx) { list(aa = g$aa[idx, ], ab = g$ab[idx, ], bb = g$bb[idx, ], n = g$n, r = g$r) } cat("all loci:", nrow(g_all$aa), " theta:", round(wc_multi(g_all)$theta, 4), "\n") cat("loci kept at MAF 0.05:", sum(keep), " theta:", round(wc_multi(subset_g(g_all, keep))$theta, 4), "\n") cat("loci discarded:", sum(!keep), " theta:", round(wc_multi(subset_g(g_all, !keep))$theta, 4), "\n") ``` The truth is 0.05 and the full set returns 0.0503. Filtering at a minor allele frequency of 0.05 throws away 2036 of 4000 loci and returns 0.0547, nine percent high, while the discarded loci on their own give 0.0309. The mechanism is not mysterious: a locus whose minor allele is rare has little variance to divide up, and the sampling noise in its frequency estimate is a larger share of what is left. Filtering removes those loci and the average shifts. The consequence is a comparability rule. An `Fst` from a filtered SNP set is not on the same scale as an `Fst` from an unfiltered one, and a threshold applied to the pooled sample is not the same as a threshold applied per population. When two studies of the same species disagree about `Fst`, the filtering step is one of the first places to look, alongside the marker type argument from the second tutorial. ## Check 4: an outlier is only an outlier under a model `Fst` outlier scans look for loci more differentiated than a neutral null and call them candidates for selection. The scan is only as good as the null, and the usual null is an island model in which every population is equally related to every other. Real sampling rarely looks like that: populations fall into regions, catchments, or lineages. Simulate two groups of three demes, with real structure among the groups and very little within them, and no selection anywhere. Then run the standard scan with an island-model null calibrated to the same overall `Fst`. ```{r outliers} set.seed(474) hier <- function(n_loci, n_ind, fst_among, fst_within, p_anc) { sh1 <- p_anc * (1 - fst_among) / fst_among sh2 <- (1 - p_anc) * (1 - fst_among) / fst_among out_aa <- out_ab <- out_bb <- matrix(0L, n_loci, 6) for (grp in 1:2) { p_grp <- rbeta(n_loci, sh1, sh2) s1 <- p_grp * (1 - fst_within) / fst_within s2 <- (1 - p_grp) * (1 - fst_within) / fst_within for (d in 1:3) { col <- (grp - 1) * 3 + d p_d <- rbeta(n_loci, s1, s2) het <- rbinom(n_loci, n_ind, 2 * p_d * (1 - p_d)) rest <- n_ind - het hom_a <- rbinom(n_loci, rest, p_d^2 / (p_d^2 + (1 - p_d)^2)) out_aa[, col] <- hom_a out_ab[, col] <- het out_bb[, col] <- rest - hom_a } } list(aa = out_aa, ab = out_ab, bb = out_bb, n = n_ind, r = 6) } p_base <- runif(2000, 0.1, 0.5) obs <- hier(2000, 30, 0.08, 0.01, p_base) obs_fit <- wc_multi(obs) cat("hierarchical data, multilocus theta:", round(obs_fit$theta, 4), "\n") null_g <- sim_loci(20000, 6, 30, obs_fit$theta, sample(p_base, 20000, TRUE)) null_per <- wc_multi(null_g)$per_locus thr <- quantile(null_per, 0.99, na.rm = TRUE) cat("island-model null, 99th percentile of per-locus theta:", round(thr, 4), "\n") cat("hierarchical loci above that threshold:", round(mean(obs_fit$per_locus > thr), 4), " nominal rate: 0.01\n") cat("null spread: sd", round(sd(null_per), 4), " hierarchical spread: sd", round(sd(obs_fit$per_locus), 4), "\n") ``` Not one locus in this simulation is under selection. The scan flags 4.8 percent of them, nearly five times the nominal one percent, because the hierarchical arrangement widens the distribution of per-locus `Fst` beyond what the island model expects (standard deviation 0.067 against 0.0445). Every flagged locus would be reported as a candidate, and a gene ontology analysis of the list would find something to say about it. ```{r fig-outliers} #| fig-width: 7 #| fig-height: 4 #| fig-cap: "Per-locus Fst under an island-model null and under a hierarchical structure with the same overall Fst, with the 99th percentile of the null marked." #| fig-alt: "Two overlapping density curves; the hierarchical curve has a heavier right tail than the island-model null, and a dashed vertical line at the null's 99th percentile cuts through it." dens_df <- data.frame( theta = c(null_per, obs_fit$per_locus), model = rep(c("island-model null", "hierarchical truth"), c(length(null_per), length(obs_fit$per_locus)))) ggplot(dens_df, aes(theta, colour = model)) + geom_density(linewidth = 0.9) + geom_vline(xintercept = thr, linetype = "dashed", colour = "grey40") + coord_cartesian(xlim = c(-0.05, 0.45)) + scale_colour_manual(values = c(`island-model null` = te_pal$sage, `hierarchical truth` = te_pal$clay), name = NULL) + labs(title = "The same overall Fst, a heavier tail", subtitle = "dashed line: 99th percentile of the island-model null", x = "per-locus Weir and Cockerham estimate", y = "density") + theme_te() ``` The fix is a null that matches the sampling: a hierarchical island model with the groups declared, or a demographic model fitted to the data first, or a scan that conditions on the observed covariance among populations. None of these is a formality. An outlier scan is a comparison against a model of the population's history, so it answers "is this locus unusual given that history", and if the history is wrong the answer is about the history. ## The honest limit The four checks constrain a population genetic analysis; they do not validate it. Two things they cannot reach are worth stating plainly. The first is that no property of allele frequencies separates selection from demography at a single locus. A bottleneck, a range expansion, or an unsampled ghost population all produce loci with extreme `Fst`, and they do it without any locus being special. The evidence that promotes a candidate from a list to a result comes from outside the scan: a repeated signal in an independent population pair, a match to an environmental gradient, a function that makes sense, an experiment. The second is that all four checks assume the marker set is a fair sample of the genome. Markers chosen because they were variable in one population, or an assay developed on one lineage, tilt every statistic in this cluster in the same direction, and no internal diagnostic can see it, because the data set does not contain the loci that were never assayed. That is the same shape of problem as the omitted trait in a selection analysis: the fix is in the design, not in the diagnostics. ## Where to go next That closes the population genetics cluster: expectations for one locus in one population, `F`-statistics for many populations, the drift and migration processes that predict `Fst`, and the checks. The natural continuation is the other place where allele frequencies change on their own, the stochastic dynamics of finite populations, where the same binomial sampling that produced drift here drives fixation probabilities and evolutionary outcomes. ## References Weir BS, Cockerham CC 1984 Evolution 38(6):1358-1370 (10.1111/j.1558-5646.1984.tb05657.x) Beaumont MA, Nichols RA 1996 Proceedings of the Royal Society B 263(1377):1619-1626 (10.1098/rspb.1996.0237) Excoffier L, Hofer T, Foll M 2009 Heredity 103(4):285-298 (10.1038/hdy.2009.74) Willing EM, Dreyer C, van Oosterhout C 2012 PLoS ONE 7(8):e42649 (10.1371/journal.pone.0042649) Lotterhos KE, Whitlock MC 2014 Molecular Ecology 23(9):2178-2192 (10.1111/mec.12725) Bhatia G, Patterson N, Sankararaman S, Price AL 2013 Genome Research 23(9):1514-1521 (10.1101/gr.154831.113) Dakin EE, Avise JC 2004 Heredity 93(5):504-509 (10.1038/sj.hdy.6800545) ## Related tutorials - [Hardy-Weinberg expectations](../hardy-weinberg-and-genotype-frequencies/) - [F-statistics and population structure](../f-statistics-and-population-structure/) - [Drift, migration and isolation by distance](../drift-migration-and-isolation-by-distance/) - [Checking a finite population model](../checking-a-finite-population-model/)