Population structure in big bluestem

0. Technical details

These exercises use the data from the assigned reading to reinforce recent lecuture materials. Everything on your end is done in R with some prepared data. We will discuss some pre-prepared complementary analyses.

You need R and nothing else.
All of the inputs are in data/ (about 0.5 MB of genotypes).
Navigate to an easy to find directory in R first
Download the data from drive

You likely need to install the following R packages. This can be done in R using the install.packages() function:

install.packages("vcfR")
install.packages("adegenet")
install.packages("poppr")
install.packages("hierfstat")
install.packages("pegas")
install.packages("vegan")
install.packages("StAMPP")
install.packages("geosphere")
install.packages("dplyr")
install.packages("ggplot2")

Two steps are shown but not run, because installing their software will not be possible for some computers. These steps include ADMIXTURE (the Q matrices are provided) and FEEMS (maps are pre-rendered). The code for both is here so you can see exactly what produced the results, but all that is required for completion is finishing the R component.

1. The dataset and the question

Andropogon gerardii (big bluestem) is a dominant prairie grass across North America. It has two common cytotypes: hexaploid (6x) and enneaploid (9x). Despite the high ploidies, both reproduce sexually and do not inbreed.

The data are from:

McAllister, C.A. & Miller, A.J. (2016) Single nucleotide polymorphism discovery via genotyping by sequencing to assess population genetic structure and recurrent polyploidization in Andropogon gerardii. American Journal of Botany 103(7): 1314–1325.

190 individuals from 24 localities, genotyped by sequencing and called against the Sorghum bicolor reference. Some quality filtering was done after acquiring the original data, which resulted in: 175 individuals from 24 localities. This includes 92 6x individuals and 83 9x individuals.

The paper’s central claim and support for that claim will be discussed Thursday. Here, we are reinforcing some basic concepts about genetic diversity and population structure. For now, let’s simply ask Is there population structure in Andropogon gerardii across North America?

A note about ploidy

We have mentioned but not discussed ploidy in depth yet. These plants are 6x and 9x, but genotypes are called as if they were diploid (0/0, 0/1, 1/1). Allele dosage is not accounted for here, and in truth, many studies ignore it. There are some approaches to deal with it directly but beyond today’s scope.

2. Where the data came from

You are starting from a cleaned Variant Call Format (VCF) file. You can find it at data/andropogon.sub.vcf.gz.

The original VCF had 397,363 single nucleotide polymorphisms (SNPs), also known as segregating sites. We reduced some of the noise by: - removing SNPs missing in more than 20% of individuals (--geno 0.2 in PLINK) - removing individuals missing more than 50% of SNPs (--mind 0.5 in PLINK) - removing near-monomorphic sites with minor allele frequency < 0.05 (--maf 0.05 in PLINK)

This left 60,141 SNPs across 175 individuals. To make some excercises run quickly, a random 10% subset of the SNPs is provided as data/andropogon.sub.vcf.gz.

Let’s explore what a VCF is together.

3. Reading the data into R

library(vcfR)
library(adegenet)
library(poppr)
library(hierfstat)
library(pegas)
library(vegan)
library(StAMPP)
library(geosphere)
library(dplyr)
library(ggplot2)

Read the VCF once and derive two objects from it. adegenet uses two representations and they are not interchangeable:

genlight — compact, one byte per genotype. Used for PCA and $F_{ST}$.
genind — allele-count table, heavier but required by poppr and hierfstat for AMOVA.

Both objects show the basic representation of population genetic data. We have counts of alleles and these can be represented as integers.

vcf <- read.vcfR("data/andropogon.sub.vcf.gz")

gl <- vcfR2genlight(vcf)

ploidy(gl) <- 2

n_individuals <- nInd(gl)
n_snps <- nLoc(gl)

cat(paste(n_individuals, "individuals and", n_snps, "SNPs\n"))

Now attach the metadata. Never assume the metadata file and the genotype file are in the same order. This is a common issue with R or other analyses. We need to relate multiple

meta <- read.csv("data/andropogon_metadata.csv", stringsAsFactors = FALSE)

# Reorder metadata to match the genotype object, then verify. If this check ever
# fails, every result below would be silently wrong -- individuals would carry
# other individuals' localities.
meta <- meta[match(indNames(gl), meta$ID), ]
stopifnot(identical(meta$ID, indNames(gl)))

pop(gl) <- as.factor(meta$Locality)

n_localities <- length(unique(meta$Locality))
n_ploidies <- length(unique(meta$Ploidy))

cat(paste("localities: ", n_localities, "   cytotypes: ", n_ploidies, "\n"))
table(meta$Locality, meta$Ploidy)

The table already tells a lot. Some localities are pure 6x, some pure 9x, many are mixed ploidy.

4. PCA

Principal components analysis on the genotype matrix. No population labels go in; any structure that appears is coming from the genotypes alone. This is often the first step to visualizing genetic structure.

# We will only retain the first 3 axes of the PCA.
pca <- glPca(gl, nf = 3, parallel = FALSE)

# glPca returns eigenvalues; convert to percent of total variance so the axes are interpretable.
eig_pct <- 100 * pca$eig / sum(pca$eig)
round(eig_pct[1:5], 2)

ploidy_levels <- sort(unique(meta$Ploidy))
ploidy_colors <- setNames(c("#1b9e77", "#d95f02"), ploidy_levels)  # 6x, 9x

pca_df <- data.frame(
  PC1 = pca$scores[, 1], PC2 = pca$scores[, 2], PC3 = pca$scores[, 3],
  Ploidy = factor(meta$Ploidy), Locality = meta$Locality, Longitude = meta$longitude
)

ggplot(pca_df, aes(PC1, PC2, color = Ploidy)) +
  geom_point(size = 2, alpha = 0.85) +
  scale_color_manual(values = ploidy_colors, name = "Cytotype", labels = ploidy_levels) +
  geom_hline(yintercept = 0, linetype = 2, colour = "grey60") +
  geom_vline(xintercept = 0, linetype = 2, colour = "grey60") +
  labs(x = sprintf("PC1 (%.1f%%)", eig_pct[1]),
       y = sprintf("PC2 (%.1f%%)", eig_pct[2])) +
  theme_bw()

Is a PC1 explaining only ~2% too low? I would say pretty normal for a widespread, wind-pollinated, outcrossing plant. Gene flow between population creates a scenario where most genetic variation is within populations rather than between them. For a selfing annual, PC1 should hit the double digits. A low PC1 is a result. It tells you structure is weak, which is why multiple approaches are needed to understand population genetic processes.

Q1: Is there a clear separation by cytotype in the PCA plot? Do those results have implications for the origins of 9x cytotypes?

5. Pairwise $F_{ST}$

$F_{ST}$ measures genetic variation between populations. Recall that 0 means panmixia and 1 means no shared variation. Recall that $F_{ST}$ is relative too.

# StAMPP computes the Weir & Cockerham estimator from the genlight object.
fst_st <- stamppFst(gl, nboots = 1, percent = 95, nclusters = 2)
fst_mat <- as.matrix(if (is.list(fst_st)) fst_st$Fsts else fst_st)

# StAMPP fills only the lower triangle; mirror it so we can subset by name later.
fst_mat[upper.tri(fst_mat)] <- t(fst_mat)[upper.tri(fst_mat)]

mean_fst <- mean(fst_mat, na.rm = TRUE)

cat(paste("mean pairwise Fst among localities:", mean_fst, "\n"))

You should get 0.021. That is quite low and it means that much of the genetic variation is shared among populations . Two populations 2,000 km apart share almost all of their genetic variation.

6. Heterozygosity by locality and cytotype

# Ho for an individual = the fraction of scored loci at which it is a heterozygote.
# Computed straight from the genlight dosage matrix (0/1/2 copies of the alt allele).
dosage <- as.matrix(gl)
meta$ind_Ho <- rowMeans(dosage == 1, na.rm = TRUE)

ggplot(meta, aes(reorder(Locality, longitude), ind_Ho)) +
  geom_boxplot(fill = "grey90", outlier.shape = NA) +
  geom_jitter(aes(color = factor(Ploidy)), width = 0.2, size = 1.5, alpha = 0.85) +
  scale_color_manual(values = ploidy_colors, name = "Cytotype",
                     labels = paste0(ploidy_levels, "x")) +
  labs(x = "Locality (west to east)", y = "Observed heterozygosity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 7))

Here, we are looking at the observed heterozygosity per individual, grouped by locality and colored by cytotype. This is not the observed heterozygosity per population per locus as calculated previously, but rather per individual. Is there any general trend across localities or cytotypes?

7. AMOVA — geography or ploidy?

Analysis of Molecular Variance (AMOVA) works like an ANOVA for genetic data: it partitions the total genetic variance into levels of a hierarchy. This allows us to measure how much of the total genetic variance is explained by cytotype versus locality.

# AMOVA needs a genind (allele counts), not a genlight.
gi <- vcfR2genind(vcf)
gi <- gi[indNames(gl), ]

pop(gi) <- as.factor(meta$Locality)
strata(gi) <- data.frame(Locality = meta$Locality, Ploidy = as.factor(meta$Ploidy))

# AMOVA cannot handle missing data; impute each locus to its mean, then recast to
# integer because the count-based statistics expect integers.
# Our PCA also imputed missing values, but that was handled in the background.
gi_imp <- missingno(gi, type = "mean")
mode(gi_imp@tab) <- "integer"

# within = FALSE is deliberate. With the default, poppr tries to partition
# WITHIN-individual variance and warns that this "cannot be calculated until the
# dosage is correctly estimated" for ambiguous-dosage data.
# For the polyploids, we can ignore the $F_{IS}$ since we are mostly interested in the between-population variance.
# We use permutations to assess the significance of the AMOVA results.
NPERM <- 9999
amova_res <- poppr.amova(gi_imp, ~Ploidy/Locality, method = "pegas",
                         nperm = NPERM, within = FALSE)
amova_res

Read the Phi-statistics:

Statistic	Meaning	Value
$\Phi_{CT}$	variance explained by cytotype	≈ 0.004
$\Phi_{SC}$	variance explained by locality within cytotype	≈ 0.045

Q2: What can you conclude about the effects of cytotype versus geography on population structure in Andropogon gerardii?

8. Isolation by distance

If gene flow simply declines with distance, genetic distance should increase with geographic distance. That is isolation by distance (IBD), and it is the null model the migration surface in §10 is measured against.

sites <- meta %>% distinct(Locality, longitude, latitude) %>% arrange(Locality)
# We need to attach the geographic data to our Fst matrix
common <- intersect(sites$Locality, rownames(fst_mat))
sites <- sites %>% filter(Locality %in% common) %>% arrange(match(Locality, common))

# Great-circle distance in km between every pair of localities.
geo_km <- distm(as.matrix(sites[, c("longitude", "latitude")]), fun = distGeo) / 1000

# Rousset's linearisation: Fst/(1-Fst) is what should be linear in distance.
fst_lin <- fst_mat[common, common] / (1 - fst_mat[common, common])
fst_lin[!is.finite(fst_lin)] <- NA

d_geo <- as.dist(geo_km)
d_gen <- as.dist(fst_lin)

# Mantel tests correlation between two DISTANCE matrices. An ordinary correlation
# test would be invalid: the entries are not independent (they share populations),
# so significance comes from permutation instead.
mantel_res <- mantel(d_gen, d_geo, permutations = NPERM, na.rm = TRUE)
mantel_res

ggplot(data.frame(distance_km = as.vector(d_geo), fst_lin = as.vector(d_gen)),
       aes(distance_km, fst_lin)) +
  geom_point(alpha = 0.4, size = 1.2) +
  geom_smooth(method = "lm", se = TRUE, colour = "#d95f02") +
  labs(x = "Geographic distance (km)", y = expression(F[ST] / (1 - F[ST])),
       title = sprintf("Mantel r = %.3f, p = %.4f",
                       mantel_res$statistic, mantel_res$signif)) +
  theme_bw()

You should get r ≈ 0.49, p = 0.001.

Q3: Does the Mantel test support isolation by distance in Andropogon gerardii?

IBD makes some simplifying assumptions that the landscape is uniform and only distance matters. Some following methods can help address departures from those assumptions, but we will do them together.

9. ADMIXTURE — where the Q matrix comes from

Ancestry-proportion analysis assigns each individual a set of proportions across $K$ ancestral components. Those proportions form the Q matrix: one row per individual, $K$ columns, each row summing to 1.

The Q matrices are provided in data/admixture/. This is the command that made them, but it requires installing the admixutre software:

# K = 1..8, with 5-fold cross-validation to choose K.
for K in $(seq 1 8); do
  admixture --cv=5 -j2 andropogon.admix.bed $K
done

Choosing K

cv <- read.table("data/admixture/cv_error.txt", col.names = c("K", "cv_error"))
best_k <- cv$K[which.min(cv$cv_error)]

ggplot(cv, aes(K, cv_error)) +
  geom_line() +
  geom_point(size = 2) +
  geom_point(data = cv[cv$K == best_k, ], size = 5, shape = 21,
             colour = "#d95f02", stroke = 1.5) +
  scale_x_continuous(breaks = cv$K) +
  labs(x = "K (number of ancestral components)", y = "cross-validation error",
       title = sprintf("Lowest CV error at K = %d", best_k)) +
  theme_bw()

The curve barely dips. With $F_{ST}$ = 0.02 there is not much structure for ADMIXTURE to find, so no K is dramatically better than its neighbours.

Some cases show a sharp optimum. Real data often looks like this, and reporting ambiguity rather than picking one K that looks the best is the best approach.

The Q matrix as a picture

PALETTE <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a",
             "#66a61e", "#e6ab02", "#a6761d", "#666666")

# Order individuals by longitude so the plot reads west (left) to east (right),
# matching how the map in section 10 is laid out.
ord <- order(meta$longitude)

# Reshape every Q matrix into one long data frame: individual x component x K.
q_long <- do.call(rbind, lapply(2:8, function(k) {
  Q <- as.matrix(read.table(sprintf("data/admixture/andropogon.admix.%d.Q", k)))
  Q <- Q[ord, , drop = FALSE]
  data.frame(
    ind       = rep(seq_len(nrow(Q)), times = k),
    component = factor(rep(seq_len(k), each = nrow(Q))),
    prop      = as.vector(Q),
    K         = factor(sprintf("K = %d%s", k, ifelse(k == best_k, " *", "")),
                       levels = sprintf("K = %d%s", 2:8, ifelse(2:8 == best_k, " *", "")))
  )
}))

# Locality boundaries, for vertical separators
loc <- meta$Locality[ord]
bounds <- which(loc[-1] != loc[-length(loc)]) + 0.5
centres <- (c(0, bounds) + c(bounds, length(loc))) / 2
labels  <- loc[round(centres)]

ggplot(q_long, aes(ind, prop, fill = component)) +
  geom_col(width = 1) +
  facet_wrap(~K, ncol = 1, strip.position = "left") +
  geom_vline(xintercept = bounds, colour = "white", linewidth = 0.3) +
  scale_fill_manual(values = PALETTE, guide = "none") +
  scale_x_continuous(breaks = centres, labels = labels, expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(x = NULL, y = NULL,
       title = "Individuals ordered west to east;  * = best K by cross-validation") +
  theme_minimal(base_size = 9) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6),
        axis.text.y = element_blank(), panel.grid = element_blank(),
        strip.placement = "outside", strip.text.y.left = element_text(angle = 0))

At K = 2 the split is west versus east: the New Mexico and Colorado sites carry one component, everything from Kansas eastward carries the other. Higher K progressively fragments that pattern without a clear geographic story.

10. FEEMS — the migration surface

FEEMS is a Python package with a fragile dependency stack (it needs NumPy < 2 and scikit-sparse < 0.5, among others). We are not installing it in class. The code below is real and runnable once FEEMS is installed.

Everything so far has been aspatial or, at best, one-dimensional (distance). FEEMS (Fast Estimation of Effective Migration Surfaces) fits a grid over the landscape and estimates, for every edge in that grid, how easily genes move across it. The result is a map:

Orange = migration lower than expected under plain IBD — a barrier
Blue = migration higher than expected — a corridor

# what produces the map (see 04_run_feems.py)
sp_graph = SpatialGraph(genotypes, coord, grid, edges, scale_snps=True)
sp_graph.fit(lamb=2.0, lamb_q=10.0, optimize_q="n-dim")

v = Viz(ax, sp_graph, projection=projection, edge_width=0.5)
v.draw_map()
v.draw_edges(use_weights=True)
v.draw_obs_nodes(use_ids=False)

The fitted migration surface. 175 individuals, 60,141 SNPs, 23 sampled demes.

There is a clear north–south corridor of reduced gene flow through the western Great Plains, and an extensive high-migration region to the east.

For discussion Thursday: Do these results support or contradict the authors’ conclusions?

The Q matrix, drawn on the map

Handing FEEMS the same ADMIXTURE Q matrix you plotted for K=2 draws it as a pie chart at each deme instead of a bar at each individual:

The same surface with ADMIXTURE K = 2 ancestry as pie charts. Compare with the K = 2 row of the barplot in section 9 — identical numbers, different geometry.

The west/east ancestry split you saw as coloured bars is now visibly sitting on top of the barrier FEEMS inferred independently from the genotypes. Two different methods, same boundary.

Cytotype on the same surface

We can draw any per-individual proportion matrix this way. Making a two-column matrix of “is this plant 6x or 9x” gives the cytotype composition of every deme:

Cytotype composition per deme (green = 6x, orange = 9x) on the migration surface.

Now compare this with the ancestry map above. The cytotypes are mosaic, not split by the barrier.

Q4 How does this agree with your AMOVA analyses?

11. Takeaways

Weak structure is a finding. PC1 ≈ 2%, $F_{ST}$ ≈ 0.02: big bluestem is a single, well-mixed, continental population.
Geography beats ploidy, by about tenfold. $\Phi_{SC}$ ≈ 0.045 vs $\Phi_{CT}$ ≈ 0.004. The 9x cytotype does not form a genetic group, which is what recurrent origins predict.
IBD is real but incomplete. Mantel r ≈ 0.49. Distance matters, but the landscape is not uniform, and FEEMS shows where.
A Q matrix is just a table of proportions. Expressed pre-individual or per-population, you can visualize ancestry proportions.

Questions to discuss

The low-migration corridor runs north–south through the western plains. What might cause it — and how would you tell a biological barrier from an artefact of where people happened to collect samples?
Does the PCA give you any information about geographic boundaries?
Cross-validation picked K = 2, but the curve is nearly flat. What would you report in a paper, and what would you need to distinguish K = 2 from K = 3?
The Mantel test says distance predicts differentiation. Does that make the FEEMS surface redundant? What does it add?

--- title: "Population structure in big bluestem" format: html: toc: true toc-depth: 2 theme: cosmo code-fold: false code-tools: true embed-resources: true fig-width: 7 fig-height: 5 --- # 0. Technical details These exercises use the data from the assigned reading to reinforce recent lecuture materials. Everything on your end is done in R with some prepared data. We will discuss some pre-prepared complementary analyses. - **You need R and nothing else.** - All of the inputs are in `data/` (about 0.5 MB of genotypes). - Navigate to an easy to find directory in R first - Download the data from [drive](https://drive.google.com/drive/folders/1zsFKlWaN5SXFxYguKL_AAeQZUc4XYasA?usp=drive_link) You likely need to install the following R packages. This can be done in R using the `install.packages()` function: ```{r} #| label: packages #| eval: false #| message: false #| warning: false install.packages("vcfR") install.packages("adegenet") install.packages("poppr") install.packages("hierfstat") install.packages("pegas") install.packages("vegan") install.packages("StAMPP") install.packages("geosphere") install.packages("dplyr") install.packages("ggplot2") ``` Two steps are **shown but not run**, because installing their software will not be possible for some computers. These steps include ADMIXTURE (the Q matrices are provided) and FEEMS (maps are pre-rendered). The code for both is here so you can see exactly what produced the results, but all that is required for completion is finishing the R component. # 1. The dataset and the question *Andropogon gerardii* (big bluestem) is a dominant prairie grass across North America. It has two common cytotypes: **hexaploid (6x)** and **enneaploid (9x)**. Despite the high ploidies, both reproduce sexually and do not inbreed. The data are from: > McAllister, C.A. & Miller, A.J. (2016) Single nucleotide polymorphism discovery via genotyping by > sequencing to assess population genetic structure and recurrent polyploidization in > *Andropogon gerardii*. *American Journal of Botany* 103(7): 1314–1325. 190 individuals from 24 localities, genotyped by sequencing and called against the *Sorghum bicolor* reference. Some quality filtering was done after acquiring the original data, which resulted in: **175 individuals from 24 localities**. This includes 92 6x individuals and 83 9x individuals. The paper's central claim and support for that claim will be discussed Thursday. Here, we are reinforcing some basic concepts about genetic diversity and population structure. For now, let's simply ask **Is there population structure in *Andropogon gerardii* across North America?** ::: {.callout-important} ## A note about ploidy We have mentioned but not discussed ploidy in depth yet. These plants are 6x and 9x, but genotypes are called as if they were **diploid** (`0/0`, `0/1`, `1/1`). Allele *dosage* is not accounted for here, and in truth, many studies ignore it. There are some approaches to deal with it directly but beyond today's scope. ::: # 2. Where the data came from You are starting from a cleaned **Variant Call Format (VCF)** file. You can find it at `data/andropogon.sub.vcf.gz`. The original VCF had 397,363 single nucleotide polymorphisms (SNPs), also known as segregating sites. We reduced some of the noise by: - removing SNPs missing in more than 20% of individuals (`--geno 0.2` in PLINK) - removing individuals missing more than 50% of SNPs (`--mind 0.5` in PLINK) - removing near-monomorphic sites with minor allele frequency < 0.05 (`--maf 0.05` in PLINK) This left **60,141 SNPs across 175 individuals**. To make some excercises run quickly, a random **10%** subset of the SNPs is provided as `data/andropogon.sub.vcf.gz`. Let's explore what a VCF is together. # 3. Reading the data into R ```{r} #| label: setup #| eval: false #| message: false #| warning: false library(vcfR) library(adegenet) library(poppr) library(hierfstat) library(pegas) library(vegan) library(StAMPP) library(geosphere) library(dplyr) library(ggplot2) ``` Read the VCF once and derive two objects from it. `adegenet` uses two representations and they are not interchangeable: - **`genlight`** — compact, one byte per genotype. Used for PCA and $F_{ST}$. - **`genind`** — allele-count table, heavier but required by `poppr` and `hierfstat` for AMOVA. Both objects show the basic representation of population genetic data. We have counts of alleles and these can be represented as integers. ```{r} #| label: read-vcf #| eval: false #| message: false #| warning: false vcf <- read.vcfR("data/andropogon.sub.vcf.gz") gl <- vcfR2genlight(vcf) ploidy(gl) <- 2 n_individuals <- nInd(gl) n_snps <- nLoc(gl) cat(paste(n_individuals, "individuals and", n_snps, "SNPs\n")) ``` Now attach the metadata. **Never assume the metadata file and the genotype file are in the same order**. This is a common issue with R or other analyses. We need to relate multiple ```{r} #| label: metadata #| eval: false meta <- read.csv("data/andropogon_metadata.csv", stringsAsFactors = FALSE) # Reorder metadata to match the genotype object, then verify. If this check ever # fails, every result below would be silently wrong -- individuals would carry # other individuals' localities. meta <- meta[match(indNames(gl), meta$ID), ] stopifnot(identical(meta$ID, indNames(gl))) pop(gl) <- as.factor(meta$Locality) n_localities <- length(unique(meta$Locality)) n_ploidies <- length(unique(meta$Ploidy)) cat(paste("localities: ", n_localities, " cytotypes: ", n_ploidies, "\n")) table(meta$Locality, meta$Ploidy) ``` The table already tells a lot. Some localities are pure 6x, some pure 9x, many are mixed ploidy. # 4. PCA Principal components analysis on the genotype matrix. No population labels go in; any structure that appears is coming from the genotypes alone. This is often the first step to visualizing genetic structure. ```{r} #| label: pca #| message: false #| eval: false # We will only retain the first 3 axes of the PCA. pca <- glPca(gl, nf = 3, parallel = FALSE) # glPca returns eigenvalues; convert to percent of total variance so the axes are interpretable. eig_pct <- 100 * pca$eig / sum(pca$eig) round(eig_pct[1:5], 2) ``` ```{r} #| label: pca-plot #| fig-cap: "PC1 vs PC2, coloured by cytotype. If ploidy were the main axis of variation, the two colours would separate." #| eval: false ploidy_levels <- sort(unique(meta$Ploidy)) ploidy_colors <- setNames(c("#1b9e77", "#d95f02"), ploidy_levels) # 6x, 9x pca_df <- data.frame( PC1 = pca$scores[, 1], PC2 = pca$scores[, 2], PC3 = pca$scores[, 3], Ploidy = factor(meta$Ploidy), Locality = meta$Locality, Longitude = meta$longitude ) ggplot(pca_df, aes(PC1, PC2, color = Ploidy)) + geom_point(size = 2, alpha = 0.85) + scale_color_manual(values = ploidy_colors, name = "Cytotype", labels = ploidy_levels) + geom_hline(yintercept = 0, linetype = 2, colour = "grey60") + geom_vline(xintercept = 0, linetype = 2, colour = "grey60") + labs(x = sprintf("PC1 (%.1f%%)", eig_pct[1]), y = sprintf("PC2 (%.1f%%)", eig_pct[2])) + theme_bw() ``` **Is a PC1 explaining only ~2% too low?** I would say pretty normal for a widespread, wind-pollinated, **outcrossing** plant. Gene flow between population creates a scenario where most genetic variation is *within* populations rather than between them. For a selfing annual, PC1 should hit the double digits. A low PC1 is a result. It tells you structure is weak, which is why multiple approaches are needed to understand population genetic processes. **Q1: Is there a clear separation by cytotype in the PCA plot? Do those results have implications for the origins of 9x cytotypes?** # 5. Pairwise $F_{ST}$ $F_{ST}$ measures genetic variation *between* populations. Recall that 0 means panmixia and 1 means no shared variation. Recall that $F_{ST}$ is relative too. ```{r} #| label: fst #| message: false #| warning: false #| eval: false # StAMPP computes the Weir & Cockerham estimator from the genlight object. fst_st <- stamppFst(gl, nboots = 1, percent = 95, nclusters = 2) fst_mat <- as.matrix(if (is.list(fst_st)) fst_st$Fsts else fst_st) # StAMPP fills only the lower triangle; mirror it so we can subset by name later. fst_mat[upper.tri(fst_mat)] <- t(fst_mat)[upper.tri(fst_mat)] mean_fst <- mean(fst_mat, na.rm = TRUE) cat(paste("mean pairwise Fst among localities:", mean_fst, "\n")) ``` You should get **0.021**. That is quite low and it means that much of the genetic variation is shared among populations . Two populations 2,000 km apart share almost all of their genetic variation. # 6. Heterozygosity by locality and cytotype ```{r} #| label: het #| fig-cap: "Observed heterozygosity per individual, by locality ordered west to east." #| eval: false # Ho for an individual = the fraction of scored loci at which it is a heterozygote. # Computed straight from the genlight dosage matrix (0/1/2 copies of the alt allele). dosage <- as.matrix(gl) meta$ind_Ho <- rowMeans(dosage == 1, na.rm = TRUE) ggplot(meta, aes(reorder(Locality, longitude), ind_Ho)) + geom_boxplot(fill = "grey90", outlier.shape = NA) + geom_jitter(aes(color = factor(Ploidy)), width = 0.2, size = 1.5, alpha = 0.85) + scale_color_manual(values = ploidy_colors, name = "Cytotype", labels = paste0(ploidy_levels, "x")) + labs(x = "Locality (west to east)", y = "Observed heterozygosity") + theme_bw() + theme(axis.text.x = element_text(angle = 60, hjust = 1, size = 7)) ``` Here, we are looking at the observed heterozygosity per individual, grouped by locality and colored by cytotype. This is not the observed heterozygosity per population per locus as calculated previously, but rather per individual. Is there any general trend across localities or cytotypes? # 7. AMOVA — geography or ploidy? Analysis of Molecular Variance (AMOVA) works like an ANOVA for genetic data: it partitions the total genetic variance into levels of a hierarchy. This allows us to measure how much of the total genetic variance is explained by **cytotype** versus **locality**. ```{r} #| label: genind #| message: false #| warning: false #| eval: false # AMOVA needs a genind (allele counts), not a genlight. gi <- vcfR2genind(vcf) gi <- gi[indNames(gl), ] pop(gi) <- as.factor(meta$Locality) strata(gi) <- data.frame(Locality = meta$Locality, Ploidy = as.factor(meta$Ploidy)) # AMOVA cannot handle missing data; impute each locus to its mean, then recast to # integer because the count-based statistics expect integers. # Our PCA also imputed missing values, but that was handled in the background. gi_imp <- missingno(gi, type = "mean") mode(gi_imp@tab) <- "integer" ``` ```{r} #| label: amova #| message: false #| warning: false #| eval: false # within = FALSE is deliberate. With the default, poppr tries to partition # WITHIN-individual variance and warns that this "cannot be calculated until the # dosage is correctly estimated" for ambiguous-dosage data. # For the polyploids, we can ignore the $F_{IS}$ since we are mostly interested in the between-population variance. # We use permutations to assess the significance of the AMOVA results. NPERM <- 9999 amova_res <- poppr.amova(gi_imp, ~Ploidy/Locality, method = "pegas", nperm = NPERM, within = FALSE) amova_res ``` Read the **Phi-statistics**: | Statistic | Meaning | Value | |---|---|---| | $\Phi_{CT}$ | variance explained by **cytotype** | ≈ 0.004 | | $\Phi_{SC}$ | variance explained by **locality within cytotype** | ≈ 0.045 | **Q2: What can you conclude about the effects of cytotype versus geography on population structure in *Andropogon gerardii*?** # 8. Isolation by distance If gene flow simply declines with distance, genetic distance should increase with geographic distance. That is **isolation by distance (IBD)**, and it is the null model the migration surface in §10 is measured against. ```{r} #| label: ibd #| fig-cap: "Linearised FST against great-circle distance between localities." #| eval: false sites <- meta %>% distinct(Locality, longitude, latitude) %>% arrange(Locality) # We need to attach the geographic data to our Fst matrix common <- intersect(sites$Locality, rownames(fst_mat)) sites <- sites %>% filter(Locality %in% common) %>% arrange(match(Locality, common)) # Great-circle distance in km between every pair of localities. geo_km <- distm(as.matrix(sites[, c("longitude", "latitude")]), fun = distGeo) / 1000 # Rousset's linearisation: Fst/(1-Fst) is what should be linear in distance. fst_lin <- fst_mat[common, common] / (1 - fst_mat[common, common]) fst_lin[!is.finite(fst_lin)] <- NA d_geo <- as.dist(geo_km) d_gen <- as.dist(fst_lin) # Mantel tests correlation between two DISTANCE matrices. An ordinary correlation # test would be invalid: the entries are not independent (they share populations), # so significance comes from permutation instead. mantel_res <- mantel(d_gen, d_geo, permutations = NPERM, na.rm = TRUE) mantel_res ``` ```{r} #| label: ibd-plot #| message: false #| fig-cap: "Each point is a pair of localities." #| eval: false ggplot(data.frame(distance_km = as.vector(d_geo), fst_lin = as.vector(d_gen)), aes(distance_km, fst_lin)) + geom_point(alpha = 0.4, size = 1.2) + geom_smooth(method = "lm", se = TRUE, colour = "#d95f02") + labs(x = "Geographic distance (km)", y = expression(F[ST] / (1 - F[ST])), title = sprintf("Mantel r = %.3f, p = %.4f", mantel_res$statistic, mantel_res$signif)) + theme_bw() ``` You should get **r ≈ 0.49, p = 0.001**. **Q3: Does the Mantel test support isolation by distance in *Andropogon gerardii*?** IBD makes some simplifying assumptions that the landscape is uniform and only distance matters. Some following methods can help address departures from those assumptions, but we will do them together. # 9. ADMIXTURE — where the Q matrix comes from Ancestry-proportion analysis assigns each individual a set of proportions across $K$ ancestral components. Those proportions form the **Q matrix**: one row per individual, $K$ columns, each row summing to 1. The Q matrices are **provided** in `data/admixture/`. This is the command that made them, but it requires installing the admixutre software: ```bash # K = 1..8, with 5-fold cross-validation to choose K. for K in $(seq 1 8); do admixture --cv=5 -j2 andropogon.admix.bed $K done ``` ## Choosing K ```{r} #| label: cv #| fig-cap: "Cross-validation error by K. The minimum is the best-supported number of components." #| eval: false cv <- read.table("data/admixture/cv_error.txt", col.names = c("K", "cv_error")) best_k <- cv$K[which.min(cv$cv_error)] ggplot(cv, aes(K, cv_error)) + geom_line() + geom_point(size = 2) + geom_point(data = cv[cv$K == best_k, ], size = 5, shape = 21, colour = "#d95f02", stroke = 1.5) + scale_x_continuous(breaks = cv$K) + labs(x = "K (number of ancestral components)", y = "cross-validation error", title = sprintf("Lowest CV error at K = %d", best_k)) + theme_bw() ``` The curve barely dips. With $F_{ST}$ = 0.02 there is not much structure for ADMIXTURE to find, so no K is dramatically better than its neighbours. Some cases show a sharp optimum. Real data often looks like this, and reporting ambiguity rather than picking one K that looks the best is the best approach. ## The Q matrix as a picture ```{r} #| label: barplots #| fig-width: 11 #| fig-height: 9 #| fig-cap: "Stacked ancestry barplots. Each column is an individual; each colour an ancestral component; every column sums to 1." #| eval: false PALETTE <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e", "#e6ab02", "#a6761d", "#666666") # Order individuals by longitude so the plot reads west (left) to east (right), # matching how the map in section 10 is laid out. ord <- order(meta$longitude) # Reshape every Q matrix into one long data frame: individual x component x K. q_long <- do.call(rbind, lapply(2:8, function(k) { Q <- as.matrix(read.table(sprintf("data/admixture/andropogon.admix.%d.Q", k))) Q <- Q[ord, , drop = FALSE] data.frame( ind = rep(seq_len(nrow(Q)), times = k), component = factor(rep(seq_len(k), each = nrow(Q))), prop = as.vector(Q), K = factor(sprintf("K = %d%s", k, ifelse(k == best_k, " *", "")), levels = sprintf("K = %d%s", 2:8, ifelse(2:8 == best_k, " *", ""))) ) })) # Locality boundaries, for vertical separators loc <- meta$Locality[ord] bounds <- which(loc[-1] != loc[-length(loc)]) + 0.5 centres <- (c(0, bounds) + c(bounds, length(loc))) / 2 labels <- loc[round(centres)] ggplot(q_long, aes(ind, prop, fill = component)) + geom_col(width = 1) + facet_wrap(~K, ncol = 1, strip.position = "left") + geom_vline(xintercept = bounds, colour = "white", linewidth = 0.3) + scale_fill_manual(values = PALETTE, guide = "none") + scale_x_continuous(breaks = centres, labels = labels, expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0)) + labs(x = NULL, y = NULL, title = "Individuals ordered west to east; * = best K by cross-validation") + theme_minimal(base_size = 9) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6), axis.text.y = element_blank(), panel.grid = element_blank(), strip.placement = "outside", strip.text.y.left = element_text(angle = 0)) ``` At **K = 2** the split is west versus east: the New Mexico and Colorado sites carry one component, everything from Kansas eastward carries the other. Higher K progressively fragments that pattern without a clear geographic story. # 10. FEEMS — the migration surface FEEMS is a Python package with a fragile dependency stack (it needs NumPy < 2 and `scikit-sparse` < 0.5, among others). We are not installing it in class. The code below is real and runnable once FEEMS is installed. Everything so far has been aspatial or, at best, one-dimensional (distance). **FEEMS** (Fast Estimation of Effective Migration Surfaces) fits a grid over the landscape and estimates, for every edge in that grid, how easily genes move across it. The result is a map: - **Orange** = migration lower than expected under plain IBD — a barrier - **Blue** = migration higher than expected — a corridor ```python # what produces the map (see 04_run_feems.py) sp_graph = SpatialGraph(genotypes, coord, grid, edges, scale_snps=True) sp_graph.fit(lamb=2.0, lamb_q=10.0, optimize_q="n-dim") v = Viz(ax, sp_graph, projection=projection, edge_width=0.5) v.draw_map() v.draw_edges(use_weights=True) v.draw_obs_nodes(use_ids=False) ``` ![The fitted migration surface. 175 individuals, 60,141 SNPs, 23 sampled demes.](data/figures/andropogon_surface.png) There is a clear north–south corridor of **reduced** gene flow through the western Great Plains, and an extensive high-migration region to the east. **For discussion Thursday: Do these results support or contradict the authors' conclusions?** ## The Q matrix, drawn on the map Handing FEEMS the same ADMIXTURE Q matrix you plotted for K=2 draws it as a pie chart at each deme instead of a bar at each individual: ![The same surface with ADMIXTURE K = 2 ancestry as pie charts. Compare with the K = 2 row of the barplot in section 9 — identical numbers, different geometry.](data/figures/andropogon_K2_cytotype.png) The west/east ancestry split you saw as coloured bars is now visibly *sitting on top of* the barrier FEEMS inferred independently from the genotypes. Two different methods, same boundary. ## Cytotype on the same surface We can draw any per-individual proportion matrix this way. Making a two-column matrix of "is this plant 6x or 9x" gives the cytotype composition of every deme: ![Cytotype composition per deme (green = 6x, orange = 9x) on the migration surface.](data/figures/andropogon_cytotype.png) Now compare this with the ancestry map above. The cytotypes are **mosaic**, not split by the barrier. **Q4 How does this agree with your AMOVA analyses?** # 11. Takeaways 1. **Weak structure is a finding.** PC1 ≈ 2%, $F_{ST}$ ≈ 0.02: big bluestem is a single, well-mixed, continental population. 2. **Geography beats ploidy, by about tenfold.** $\Phi_{SC}$ ≈ 0.045 vs $\Phi_{CT}$ ≈ 0.004. The 9x cytotype does not form a genetic group, which is what recurrent origins predict. 3. **IBD is real but incomplete.** Mantel r ≈ 0.49. Distance matters, but the landscape is not uniform, and FEEMS shows where. 4. **A Q matrix is just a table of proportions.** Expressed pre-individual or per-population, you can visualize ancestry proportions. ## Questions to discuss 1. The low-migration corridor runs north–south through the western plains. What might cause it — and how would you tell a *biological* barrier from an artefact of where people happened to collect samples? 2. Does the PCA give you any information about geographic boundaries? 3. Cross-validation picked K = 2, but the curve is nearly flat. What would you report in a paper, and what would you need to distinguish K = 2 from K = 3? 4. The Mantel test says distance predicts differentiation. Does that make the FEEMS surface redundant? What does it add?