R: Population Genetics Foundations

hpc

foundations

github

popgen

Allele frequencies, heterozygosity, and a first genotype–environment association

Author

Anoob Prakash

Published

October 7, 2026

PART I

Goal

Apply your R skills to population-genetics concepts: allele frequency, heterozygosity, genetic diversity, and a first genotype/phenotype–environment association (GEA), using the genetic-parameter columns already present in the red-spruce dataset (PC1, PC2, Family_Homozygosity, Population_Homozygosity, Genetic_Diversity, Genetic_Load).

We don’t have a raw genotype (VCF) file for this dataset, so this module teaches the underlying math with a small simulated SNP matrix, then applies the same reasoning to the real, pre-computed genetic-parameter columns already in your data.

Prerequisites

1. Data required data from github

Download data

Download the raw GitHub file directly into data/raw/

Copy the code below and run it inside the cluster-exercise directory

dir.create("data/raw", recursive = TRUE, showWarnings = FALSE)
download.file(
"https://raw.githubusercontent.com/anoobvinu07/Genomic_assisted_selection/master/data/FitnessTraits_GeneticParameters_RedSpruce.txt",
      destfile = "data/raw/red_spruce_fitness_traits.txt")

The -L flag tells curl to follow redirects, and -o specifies the local output file. Verify the download before doing any analysis:

Recreate the rest of the project layout if you don’t already have it:

dir.create("results", showWarnings = FALSE)
dir.create("src", showWarnings = FALSE)

Project directory setup

This excercise assumes the project folder follows the following structure:

cluster-exercise/
 ├── data/
 │   └── raw/       # downloaded, unchanged source data
 ├── logs/          # Slurm standard output and error logs
 ├── results/       # analysis products created by R
 └── src/           # R and Slurm scripts

2. Install the packages used here:

  install.packages(c("vegan", "here"))
  # Optional, for working with real genotype/VCF files beyond this module:
  # install.packages("BiocManager")
  # BiocManager::install(c("vcfR", "adegenet"))

Learning objectives

By the end of this module you should be able to:

Explain what allele frequency, observed heterozygosity, and expected heterozygosity mean
Represent genotypes as a simple numeric matrix (0/1/2 coding) and compute allele frequencies from it
Compute observed vs. expected (Hardy-Weinberg) heterozygosity per locus and interpret the difference
Run a PCA on genetic markers with prcomp() and interpret PC1/PC2 in a population-genetics context
Test a simple genotype/phenotype–environment association using correlation and linear models
Recognize the step from a simple correlation to a multivariate method (RDA) used in real genotype-environment association work
Know when to move from a plain R matrix to dedicated packages (vcfR, adegenet) for real VCF-scale data

R for bioinformatics quick reference

Function/package	Purpose	Example
`matrix()`	Build a genotype matrix	`matrix(..., nrow, ncol)`
`sample()`	Randomly draw values, optionally with replacement and specified probabilities	`sample(0:2, 300, replace = TRUE, prob = c(0.4, 0.4, 0.2))`
`rownames`/ `colnames`	Get or assign labels for matrix rows and columns	`rownames(geno) <- paste0("ind_", 1:n_ind)`
`paste0`	Convert inputs to text and concatenate them with no separator	`paste0("ind",1:30)`
`colMeans()`	Column-wise means (used for allele freq.)	`colMeans(geno) / 2`
`table()`	Count occurrences of each value	`table(pop_1$genotype)`
`prop.table()`	Convert counts to proportions	`prop.table(table(pop_1$genotype))`
`rbind()`	Stack data frames row-wise	`rbind(pop_1, pop_2)`
`prcomp()`	Principal component analysis	`prcomp(geno, scale. = TRUE)`
`cor.test()`	Test correlation between two variables	`cor.test(x, y)`
`lm()`	Fit a linear model	`lm(y ~ x1 + x2, data = df)`
`vegan::rda()`	Redundancy analysis (multivariate GEA)	`rda(Y ~ x1 + x2, data = df)`
`vcfR::read.vcfR()`	Read a real VCF file (beyond this module)	`read.vcfR("file.vcf")`
`adegenet::genind`	Genotype container object (beyond this module)	`df2genind(...)`

Basics of population genetics via simulation

A minimal SNP genotype matrix in R

Genotypes at a biallelic SNP are usually coded as the count of the alternate allele: 0 (homozygous reference), 1 (heterozygous), 2 (homozygous alternate). Let’s simulate a toy dataset — 30 individuals, 10 loci:

set.seed(123)
n_ind  <- 30 # number of individuals
n_loci <- 10 # number of SNP loci

This code creates a simulated genotype matrix where rows are individuals, columns are SNP loci, and each cell contains a diploid genotype coded as 0, 1, or 2.

geno <- matrix(
  sample(0:2,                # sample values from 0, 1 or 2
  n_ind * n_loci,            # no. of draws 30 x 10 or size of the matrix
  replace = TRUE,            # each genotype category can be drawn repeatedly
  prob = c(0.4, 0.4, 0.2)),  # sampling probability for genotypes 0,1 and 2
  nrow = n_ind,              # number of rows - filled with ind IDs
  ncol = n_loci              # number of columns - filled with SNP IDs
)

Logic check

For a biallelic SNP, the number represents the count of the alternate allele carried by an individual:

Value	Genotypes	Biological meaning
0	0/0	Homozygous reference; zero alternate alleles
1	0/1	Heterozygous; one alternate allele
2	1/1	Homozygous alternate; two alternate alleles

So, for example, a 1 in row ind_7 and column locus_3 means individual 7 is heterozygous at SNP locus 3.

Set the probabilty of each genotype for the simlutaion run with prob argument of sample() function.

Genotype value	Sampling probability
0	0.40
1	0.40
2	0.20

rownames(geno) <- paste0("ind_", 1:n_ind)
colnames(geno) <- paste0("locus_", 1:n_loci)


geno[1:5, 1:5]

      locus_1 locus_2 locus_3 locus_4 locus_5
ind_1       1       2       0       1       0
ind_2       0       2       1       0       1
ind_3       0       0       1       1       1
ind_4       2       0       1       0       1
ind_5       2       1       2       1       1

Allele frequencies

Allele frequency is the proportion of all gene copies at a locus that are a particular allele. \[p + q =1\]

Each genotype value is the count of the alternate allele out of 2, so the mean genotype divided by 2 is the alternate allele frequency for that locus:

alt_freq <- colMeans(geno) / 2
ref_freq <- 1 - alt_freq


data.frame(locus = names(alt_freq),       # col 1 filled with locus names
           alt_freq = round(alt_freq, 3), # col 2 filled with alt freq
           ref_freq = round(ref_freq, 3)) # col 3 is filled with ref freq

            locus alt_freq ref_freq
locus_1   locus_1    0.400    0.600
locus_2   locus_2    0.383    0.617
locus_3   locus_3    0.383    0.617
locus_4   locus_4    0.367    0.633
locus_5   locus_5    0.400    0.600
locus_6   locus_6    0.317    0.683
locus_7   locus_7    0.417    0.583
locus_8   locus_8    0.400    0.600
locus_9   locus_9    0.367    0.633
locus_10 locus_10    0.450    0.550

Worked example: allele frequency by hand

Suppose one locus has 20 genotyped individuals: 12 are 0/0, 5 are 0/1, and 3 are 1/1.

Each individual carries 2 allele copies, so the denominator is $2N = 40$ gene copies:

Alternate allele copies: $(0 \times 12) + (1 \times 5) + (2 \times 3) = 11$
Reference allele copies: $(2 \times 12) + (1 \times 5) + (0 \times 3) = 29$

\[p_{alt} = \frac{11}{40} = 0.275 \qquad p_{ref} = \frac{29}{40} = 0.725 \qquad p_{alt} + p_{ref} = 1\]

The shortcut mean(geno_column) / 2 works for exactly the same reason: the mean of a 0/1/2 column is the average number of alternate alleles per individual, so dividing by 2 converts it to a per-gene-copy frequency.

Worked example: why divide by 2?

The mean of a 0/1/2 genotype column counts alleles per individual, but allele frequency counts alleles per gene copy. A diploid individual has 2 copies, hence the division by 2. Sanity checks that always hold:

alt_freq must lie between 0 and 1
alt_freq + ref_freq must equal exactly 1
a column that is entirely 0 or entirely 2 is monomorphic (no variation) and carries no information for a GEA — real pipelines filter these out along with low-frequency variants

PART II

Goals

In this exercise, you will simulate genotypes for two populations at one biallelic locus, calculate observed allele and genotype frequencies, compare the populations, and combine them into a single dataset.

Learning objectives

By the end of this section of the module, you should be able to:

Use sample() to simulate diploid genotypes (AA, Aa, and aa).
Calculate genotype frequencies from observed genotype counts.
Calculate allele frequencies (p) and (q) from genotype counts.
Compare populations with different allele frequencies and population sizes.
Combine population-level data with rbind().
Explain why combining structured populations can change apparent genotype frequencies, even when each population was generated under Hardy–Weinberg expectations.

Simulating two populations at one locus

Lets start off by creating two populations with of same size (n=20), but different allele A frequency. This exercise simulates genotypes at one biallelic locus in two populations. You will use sample() to assign genotypes, calculate genotype frequencies with table() and prop.table(), then combine populations with rbind().

For this example, the two alleles are A and a, with three possible diploid genotypes:

Genotype	Meaning
AA	Two copies of allele A
Aa	One copy of A and one copy of a
aa	Two copies of allele a

At a biallelic locus:

\[p + q = 1 \] where,
p is the frequency of allele A and;
q is the frequency of allele a.

Under Hardy–Weinberg expectations, genotype probabilities are:

\[AA = p^2\] \[Aa = 2pq\]

\[aa = q^2\]

Concept review : H-W equilibrium

What is Hardy-Weinberg equilibrium?

The key expectataion under Hardy-Weinberg equilibrium is $p^2$ for AA, $2pq$ for Aa, and $q^2$ for aa, where $p + q = 1$. In finite populations, observed values vary around thse expectations because of random sampling, variation becomes more obvious at small sample sizes.

The code below uses those probabilities to generate individuals in each population. Because genotypes are randomly sampled, observed frequencies will usually be close to but not identical to the expected frequencies.

set.seed(42)


# Population 1: 20 individuals; allele A frequency p = 0.80
p_1 <- 0.80  
q_1 <- 1 - p_1


pop_1 <- data.frame(population = "Pop_1",
                    # sample function for choosing genotypes
                    genotype = sample(c("AA", "Aa", "aa"), # possible genotype values
                      size = 20,                           # no. of individuals to simulate
                      replace = TRUE,                      # genotypes can  occur more than once
                      # p2 + 2pq + q2
                      prob = c(p_1^2, 2 * p_1 * q_1, q_1^2)) # probability of each genotype
                      )


head(pop_1) # peek into the pop_1 dataframe

  population genotype
1      Pop_1       Aa
2      Pop_1       Aa
3      Pop_1       AA
4      Pop_1       Aa
5      Pop_1       Aa
6      Pop_1       AA

# Population 2: 20 individuals; allele A frequency p = 0.30
p_2 <- 0.30
q_2 <- 1 - p_2


pop_2 <- data.frame(population = "Pop_2",
                    genotype = sample(c("AA", "Aa", "aa"),
                      size = 20,
                      replace = TRUE,
                      prob = c(p_2^2, 2 * p_2 * q_2, q_2^2))
                      )


head(pop_2) # peek into the pop_2 dataframe

  population genotype
1      Pop_2       Aa
2      Pop_2       aa
3      Pop_2       AA
4      Pop_2       AA
5      Pop_2       aa
6      Pop_2       Aa

Code breakdown

sample() randomly assigns each individual one of three possible genotypes: AA, Aa, or aa.
replace = TRUE allows genotype categories to occur repeatedly.
prob = c($p^2$, $2.p.q$, $q^2$) supplies the Hardy–Weinberg genotype probabilities.

For population 1, where p = 0.80 and q = 0.20, the expected genotype frequencies are: \[AA = 0.80^2 = 0.64\] \[Aa = 2 \times 0.80 \times 0.20 = 0.32\]
\[aa = 0.20^2 = 0.04\] Therefore, Population 1 should countain mostly AA individuals.

For Population 2, where p = 0.30 and q = 0.70, the expected genotype frequencies are: \[AA = 0.30^2 = 0.09\] \[Aa = 2 \times 0.30 \times 0.70 = 0.42\]
\[aa = 0.70^2 = 0.49\] Therefore, population2 should contain mostly aa individuals.

Calculate genotype frequencies

Use table() to count each genotype and prop.table() to convert those counts into frequencies.

# Count observed genotypes in each population
geno_count_pop1 <- table(pop_1$genotype)
geno_count_pop2 <- table(pop_2$genotype)



# Convert genotype counts to genotype frequencies  
geno_prop_pop1 <- prop.table(table(pop_1$genotype))
geno_prop_pop2 <- prop.table(table(pop_2$genotype))


# Combine individuals from the two populations
pooled_population <- rbind(pop_1, pop_2)


geno_count_pooled <- table(pooled_population$genotype)
geno_prop_pooled <- prop.table(table(pooled_population$genotype))



# Build one comparison table
geno_summary <- data.frame(
  genotype = c("AA", "Aa", "aa"),
  count_pop1 = as.integer(geno_count_pop1[c("AA", "Aa", "aa")]),
  freq_pop1  = as.numeric(geno_prop_pop1[c("AA", "Aa", "aa")]),
  count_pop2 = as.integer(geno_count_pop2[c("AA", "Aa", "aa")]),
  freq_pop2  = as.numeric(geno_prop_pop2[c("AA", "Aa", "aa")]),
  count_pooled = as.integer(geno_count_pooled[c("AA", "Aa", "aa")]),
  freq_pooled  = as.numeric(geno_prop_pooled[c("AA", "Aa", "aa")])
)


geno_summary

  genotype count_pop1 freq_pop1 count_pop2 freq_pop2 count_pooled freq_pooled
1       AA          9      0.45          2       0.1           11       0.275
2       Aa         10      0.50         10       0.5           20       0.500
3       aa          1      0.05          8       0.4            9       0.225

# The genotype frequencies should sum up to 1. 
sum(geno_prop_pop1)

[1] 1

sum(geno_prop_pop2)

[1] 1

sum(geno_prop_pooled)

[1] 1

Code breakdown

rbind() or “row” bind stacks the two population data frames row by row.
table() counts each observed genotype.

Calculate allele frequencies

Genotype frequencies describe the proportion of individuals in each genotype class. Allele frequencies instead count all gene copies. Each diploid individual carries two allele copies:

Genotype	Number of A alleles	Number of a alleles
AA	2	0
Aa	1	1
aa	0	2

\[p_{A} = \frac{2n_{AA} + n_{Aa}}{2N} \]

\[q_{a} = \frac{2n_{aa} + n_{Aa}}{2N} \] Where each AA individuals contributes two A alleles, each aa individual contributes two aa alleles, and each heterozygote contributes one copy of each allele.

To calculate the allele frequency for Population 1:

# Count each genotype in Population 1
counts_1 <- geno_count_pop1


# Number of individuals
N_1 <- nrow(pop_1)


# Count A and a allele copies
n_A_1 <- 2 * counts_1["AA"] + counts_1["Aa"]
n_a_1 <- 2 * counts_1["aa"] + counts_1["Aa"]


# Convert allele counts to allele frequencies
p_A_1 <- n_A_1 / (2 * N_1)
q_a_1 <- n_a_1 / (2 * N_1)


p_A_1

 AA 
0.7

q_a_1

 aa 
0.3

p_A_1 + q_a_1

AA 
 1

Worked example: the numbers behind Population 1

With set.seed(42), Population 1 has 9 AA, 10 Aa, 1 aa. Plugging those counts into the formulas:

\[p_{A} = \frac{(2 \times 9) + 10}{40} = \frac{28}{40} = 0.70\]

\[q_{a} = \frac{(2 \times 1) + 10}{40} = \frac{12}{40} = 0.30\]

Note that $p_A + q_a = 1$ exactly — this is not a coincidence but a bookkeeping identity: every one of the $2N = 40$ gene copies is counted exactly once on either the numerator side of $p_A$ or of $q_a$.

Two things worth noticing:

The observed allele frequency (0.70) differs from the simulation input (0.80) — random sampling in a finite population (N = 20) shifts the estimate. Smaller samples drift further from the input value.
The observed genotype counts (9/10/1) do not match the H–W expectations for $p = 0.70$ (12.25/6.00/1.75 out of 20) — here we observe an excess of heterozygotes. With only 20 individuals, deviations like this are expected sampling noise, not biology.

Worked example: expected vs. observed counts

Expected counts are just $N \times$ expected frequency. For Population 1 ($p = 0.80$, $N = 20$):

Genotype	Expected frequency	Expected count (of 20)	Observed count (seed 42)
AA	0.64	12.8	9
Aa	0.32	6.4	10
aa	0.04	0.8	1

The expected aa count is less than 1 individual. This is why formal goodness-of-fit tests (like chisq.test()) are unreliable here: the standard rule of thumb requires expected counts of at least ~5 per category. With $q^2 = 0.04$ and N = 20, you would need roughly $5 / 0.04 = 125$ individuals just to expect 5 aa individuals. Rare genotype classes need large samples.

Excercises

Exercise 1 — Allele frequencies for Population 2

Estimate the allele frequencies $p_A$ and $q_a$ for Population 2, following the same steps used for Population 1: count genotypes, count allele copies, and divide by $2N$. Verify that $p + q = 1$.

Solution

counts_2 <- geno_count_pop2
N_2 <- nrow(pop_2)

n_A_2 <- 2 * counts_2["AA"] + counts_2["Aa"]
n_a_2 <- 2 * counts_2["aa"] + counts_2["Aa"]

p_A_2 <- unname(as.numeric(n_A_2 / (2 * N_2)))
q_a_2 <- unname(as.numeric(n_a_2 / (2 * N_2)))

p_A_2

[1] 0.35

q_a_2

[1] 0.65

p_A_2 + q_a_2

[1] 1

With set.seed(42), Population 2 has 2 AA, 10 Aa, 8 aa, so:

\[p_{A} = \frac{(2 \times 2) + 10}{40} = \frac{14}{40} = 0.35\]

\[q_{a} = \frac{(2 \times 8) + 10}{40} = \frac{26}{40} = 0.65\]

Again the observed $p$ (0.35) sits above the simulation input (0.30) — finite-sample noise. The two populations are strongly differentiated at this locus: $p_A$ = 0.70 vs. 0.35.

Exercise 2 — Allele frequencies for the pooled population

Compute the allele frequencies for the pooled population (Pop_1 + Pop_2). Compare them to the within-population values, and explain the pattern: where does the pooled estimate fall, and why?

Solution

counts_p <- geno_count_pooled
N_p <- nrow(pooled_population)

n_A_p <- 2 * counts_p["AA"] + counts_p["Aa"]
n_a_p <- 2 * counts_p["aa"] + counts_p["Aa"]

p_A_p <- unname(as.numeric(n_A_p / (2 * N_p)))
q_a_p <- unname(as.numeric(n_a_p / (2 * N_p)))

p_A_p

[1] 0.525

q_a_p

[1] 0.475

Pooled counts are 11 AA, 20 Aa, 9 aa (40 individuals, 80 gene copies):

\[p_{A} = \frac{(2 \times 11) + 20}{80} = \frac{42}{80} = 0.525\]

\[q_{a} = \frac{(2 \times 9) + 20}{80} = \frac{38}{80} = 0.475\]

The pooled $p_A$ = 0.525 is almost exactly the average of the two within-population frequencies:

\[\frac{p_{A,1} + p_{A,2}}{2} = \frac{0.70 + 0.35}{2} = 0.525\]

(With unequal population sizes, the pooled frequency is a weighted average, weighted by each population’s size.) Mixing does not create or destroy alleles — it just averages them. But as the next exercise shows, it does distort genotype frequencies.

Exercise 3 — Pooling changes genotype frequencies

Using the geno_summary table you built earlier:

Compute the expected genotype frequencies for the pooled population from its allele frequencies ($p^2$, $2pq$, $q^2$).
Compare them to the observed pooled genotype frequencies.
Compute the expected genotype frequencies you would get by averaging the two populations’ own H–W expectations.
Which comparison do the observed pooled frequencies resemble more? Why?

Solution

# 1. H-W expectations from the POOLED allele frequencies
exp_pooled_from_pool_p <- c(AA = p_A_p^2, Aa = 2 * p_A_p * q_a_p, aa = q_a_p^2)

# 2. Observed pooled genotype frequencies
obs_pooled <- as.numeric(geno_prop_pooled[c("AA", "Aa", "aa")])

# 3. Average of the two populations' OWN H-W expectations
exp_avg_of_pops <- c(
  AA = (p_1^2 + p_2^2) / 2,
  Aa = (2 * p_1 * q_1 + 2 * p_2 * q_2) / 2,
  aa = (q_1^2 + q_2^2) / 2
)

rbind(
  observed            = obs_pooled,
  HW_from_pooled_p    = exp_pooled_from_pool_p,
  avg_of_pop_HW       = exp_avg_of_pops
)

                       AA      Aa       aa
observed         0.275000 0.50000 0.225000
HW_from_pooled_p 0.275625 0.49875 0.225625
avg_of_pop_HW    0.365000 0.37000 0.265000

Using the seed-42 numbers ($p_{pool} = 0.525$):

H–W from pooled $p$: $AA = 0.276$, $Aa = 0.499$, $aa = 0.226$
Average of the two populations’ own expectations: $AA = 0.365$, $Aa = 0.370$, $aa = 0.265$

Here is the interesting part: at this seed, the observed pooled frequencies (0.275 / 0.500 / 0.225) land almost exactly on the H–W line from the pooled $p$ — the Wahlund deficit is invisible. Compare the two reference rows to see why the deficit should be there: H–W from the pooled $p$ predicts 49.9% heterozygotes, but heterozygotes in the mixture can only come from the two populations’ own heterozygote proportions, which average 37%. The systematic expectation is a 13-point heterozygote deficit — the Wahlund effect.

Why can’t we see it? At $N = 20$ per population, the standard error of a genotype proportion near 0.5 is about $\sqrt{0.25/20} \approx 0.11$ — the same size as the deficit itself. In this seed, Pop_1 happened to draw a heterozygote excess (10 observed vs. 6.4 expected), which cancels the Wahlund gap almost perfectly.

The takeaway is not that the effect is fake, but that at small N it is buried in noise. Pooling two H–W populations systematically produces a heterozygote deficit relative to H–W expectations computed from the pooled allele frequencies — no population evolved, only the reference frame changed. Part III re-runs this exact comparison at $N = 500$ per population, where the deficit stands out unmistakably.

Exercise 4 — Allele frequencies from the genotype matrix

Return to the 30 individuals $\times$ 10 loci matrix geno from Part I.

Compute the alternate allele frequency for every locus by hand using allele counting (2*n2 + n1 divided by 2*N), not colMeans().
Verify that your result matches colMeans(geno) / 2 exactly.
Compute the minor allele frequency (MAF) for each locus. Which locus is closest to fixation? Which is most variable?

Solution

N <- nrow(geno)

# 1. Allele counting, locus by locus
n_alt <- 2 * colSums(geno == 2) + colSums(geno == 1)
alt_freq_manual <- n_alt / (2 * N)

# 2. Cross-check with the shortcut
all.equal(as.numeric(alt_freq_manual), as.numeric(alt_freq))

[1] TRUE

# 3. Minor allele frequency
maf <- pmin(alt_freq, 1 - alt_freq)
sort(round(maf, 3))

 locus_6  locus_4  locus_9  locus_2  locus_3  locus_1  locus_5  locus_8 
   0.317    0.367    0.367    0.383    0.383    0.400    0.400    0.400 
 locus_7 locus_10 
   0.417    0.450

all.equal() confirms the two computations agree — colMeans(geno)/2 is just algebraically compressed allele counting: the mean of a 0/1/2 column is $(n_{het} + 2 \cdot n_{homAlt}) / N$, which divided by 2 is exactly the alternate allele frequency.

For the MAF, the least variable locus is locus_6 (alt = 0.317, MAF $\approx$ 0.317) and the most variable is locus_10 (alt = 0.450, MAF $\approx$ 0.450). The most informative loci for GEA are the ones with MAF near 0.5, which is why filtering pipelines typically drop loci with MAF below 0.05: rare variants carry little signal and inflate noise.

Exercise 5 — Explore the simulation

Change one thing at a time and re-run the whole Part II pipeline. Before you run each variant, write down what you expect to happen; then check.

Sample size: set both populations to size = 500. Do observed frequencies land closer to the H–W expectations?
Similar populations: set $p_2 = 0.75$. What happens to the pooled genotype frequencies and the Wahlund gap?
Divergent populations: set $p_2 = 0.05$. What happens?
Different seed: change set.seed(42) to another number. Which quantities change (counts, frequencies) and which relationships stay fixed ($p + q = 1$, pooled $p$ = average of within-population $p$)?

Solution — what you should see

Sample size (N = 500): observed genotype frequencies converge on the H–W expectations. Sampling noise shrinks like $1/\sqrt{N}$ — ten times the sample, about three times less noise. This is the single most reliable way to distinguish “sampling noise” from “real deviation from H–W”.

Similar populations ($p_2 = 0.75$): the two populations are nearly indistinguishable at this locus. Pooling barely changes anything — there is almost no Wahlund effect because there is almost no structure. Weak structure = weak heterozygote deficit.

Divergent populations ($p_2 = 0.05$): the contrast is extreme — Pop_1 is nearly all AA, Pop_2 nearly all aa. The pooled sample looks like a classic heterozygote-deficient population even though each source population is in perfect H–W equilibrium. Strong structure = strong apparent deficit. This is the situation that most often fools naive analyses of mixed samples (e.g., pooling before computing heterozygosity in a VCF).

Different seed: the counts and frequencies all shift slightly, but the relationships hold: $p + q = 1$ in every population and the pooled $p$ always equals the (size-weighted) average of within-population $p$. Identities survive randomness; estimates don’t.

One practical gotcha when you change seeds or parameters: if a genotype class happens not to occur (e.g., no aa individuals in Pop_1), table() silently drops that level and indexing counts["aa"] returns NA. Keep all three classes with table(factor(pop_1$genotype, levels = c("AA", "Aa", "aa"))) so counts of zero are explicit.

PART III

Observed and expected heterozygosity

Observed heterozygosity ($H_{o}$): the proportion of individuals with genotype 1 (heterozygous) at a locus.

Expected heterozygosity ($H_{e}$): what you’d expect under Hardy-Weinberg equilibrium, given the allele frequencies: $H_{e} = 2pq$, where p and q are the reference and alternate allele frequencies.

Concept review: What is Hardy-Weinberg equilibrium?

Ho <- colMeans(geno == 1)
He <- 2 * ref_freq * alt_freq


het_table <- data.frame(locus = names(Ho), 
                        Ho = round(Ho, 3), 
                        He = round(He, 3),
                        diff = round(Ho - He, 3))
het_table

            locus    Ho    He   diff
locus_1   locus_1 0.267 0.480 -0.213
locus_2   locus_2 0.500 0.473  0.027
locus_3   locus_3 0.367 0.473 -0.106
locus_4   locus_4 0.333 0.464 -0.131
locus_5   locus_5 0.467 0.480 -0.013
locus_6   locus_6 0.433 0.433  0.001
locus_7   locus_7 0.367 0.486 -0.119
locus_8   locus_8 0.400 0.480 -0.080
locus_9   locus_9 0.333 0.464 -0.131
locus_10 locus_10 0.567 0.495  0.072

If Ho is consistently lower than He across loci, that’s a signature of inbreeding or population substructure — the same underlying idea behind the Family_Homozygosity and Population_Homozygosity columns already computed in the real world dataset used in the previous modules.

Worked example: heterozygosity by hand

Take the hand-calculation locus from Part I: 12 0/0, 5 0/1, 3 1/1 ($N = 20$, $p_{alt} = 0.275$, $p_{ref} = 0.725$).

Observed heterozygosity — just count the heterozygotes:

\[H_{o} = \frac{n_{het}}{N} = \frac{5}{20} = 0.25\]

Expected heterozygosity — the H–W heterozygote frequency at those allele frequencies:

\[H_{e} = 2pq = 2 \times 0.725 \times 0.275 = 0.399\]

This locus has fewer heterozygotes than expected ($H_o < H_e$). The standard way to summarize the gap is the inbreeding coefficient:

\[F = 1 - \frac{H_{o}}{H_{e}} = 1 - \frac{0.25}{0.399} = 0.373\]

$F = 0$: the population matches H–W exactly at this locus
$F > 0$: heterozygote deficit (inbreeding, substructure, or chance)
$F < 0$: heterozygote excess (e.g., overdominance, or sampling noise)

Notice how the pieces from Part II assemble here: allele counts → allele frequencies → H–W expectations → observed vs. expected → F. Every downstream population-genetics quantity is built from the same two inputs, $p$ and $N$.

Worked example: $H_e$ is maximized at $p = 0.5$

$H_e = 2pq$ is a parabola in $p$ — it peaks when both alleles are equally common:

$p$	$q$	$H_e = 2pq$
0.1	0.9	0.18
0.3	0.7	0.42
0.5	0.5	0.50
0.7	0.3	0.42
0.9	0.1	0.18
1.0	0.0	0.00

Three consequences worth internalizing:

Near-fixed loci are uninformative. At $p = 0.95$, $H_e = 0.095$ — almost no variation to work with. This is why MAF filtering exists (Exercise 4, Part II).
Heterozygosity measures evenness, not richness. A population fixed for a different allele at every locus and one fixed for the same allele at every locus both have $H_e = 0$.
Homozygosity and diversity are two views of the same quantity: for a biallelic locus, $H_e + \text{homozygosity} = 1$. In the red-spruce data, this is why Population_Homozygosity and Genetic_Diversity should be near-perfect mirror images of each other (negatively correlated) — check it in the practical exercises.

Worked example: reading the simulated het_table

For the Part I simulation (seed 123), the table shows:

Mean $H_o \approx 0.403$, mean $H_e \approx 0.473$, so $H_o < H_e$ at most loci, and the overall $F \approx 0.147$.
The largest deficit is at locus_1 ($H_o = 0.267$ vs. $H_e = 0.480$).
The largest excess is at locus_10 ($H_o = 0.567$ vs. $H_e = 0.495$).

Here is the subtle part — this deficit is built into the simulation, not a biological process. We drew genotypes with probabilities c(0.4, 0.4, 0.2) rather than drawing alleles and pairing them randomly. Those genotype probabilities imply an allele frequency of $p = (0 \times 0.4 + 1 \times 0.4 + 2 \times 0.2)/2 = 0.4$, and H–W at that frequency would give $H_e = 2(0.4)(0.6) = 0.48$ — but our draw only produces heterozygotes 40% of the time. So the simulation has an inbreeding coefficient of about $F = 1 - 0.40/0.48 \approx 0.17$ by construction.

The lesson: always ask how a simulation generates its data before interpreting the output. A genotype-frequency simulator and an allele-frequency simulator with the same $p$ produce different heterozygosities — which is precisely the $H_o$ vs $H_e$ distinction.

The Wahlund effect: connecting Part II to Part III

In Part II you saw that pooling two populations distorts genotype frequencies. Here is the same demonstration at a sample size where the pattern is unambiguous. This time, simulate two large populations (N = 500 each) at one locus:

set.seed(7)
N_w <- 500

# Same structure as Part II, just bigger
pop_big_1 <- sample(c("AA", "Aa", "aa"), N_w, replace = TRUE,
                    prob = c(p_1^2, 2 * p_1 * q_1, q_1^2))   # p = 0.80
pop_big_2 <- sample(c("AA", "Aa", "aa"), N_w, replace = TRUE,
                    prob = c(p_2^2, 2 * p_2 * q_2, q_2^2))   # p = 0.30

pool_big <- c(pop_big_1, pop_big_2)

Ho_of <- function(g) sum(g == "Aa") / length(g)
p_of  <- function(g) (2 * sum(g == "AA") + sum(g == "Aa")) / (2 * length(g))

data.frame(
  group = c("Pop_1", "Pop_2", "Pooled"),
  p_A = round(c(p_of(pop_big_1), p_of(pop_big_2), p_of(pool_big)), 3),
  Ho = round(c(Ho_of(pop_big_1), Ho_of(pop_big_2), Ho_of(pool_big)), 3),
  He_from_p = round(2 * c(p_of(pop_big_1), p_of(pop_big_2), p_of(pool_big)) *
                       (1 - c(p_of(pop_big_1), p_of(pop_big_2), p_of(pool_big))), 3)
)

   group   p_A    Ho He_from_p
1  Pop_1 0.800 0.332     0.320
2  Pop_2 0.307 0.418     0.426
3 Pooled 0.553 0.375     0.494

Reading the Wahlund table

With seed 7 you should see numbers like these:

Group	$p_A$	$H_o$	$H_e$ (from its own $p$)
Pop_1	0.800	0.332	0.320
Pop_2	0.307	0.418	0.426
Pooled	0.553	0.375	0.494

Read the last row carefully — both comparisons matter:

Within each population, $H_o \approx H_e$. Each population is in H–W equilibrium, as designed. No inbreeding, no structure — nothing is wrong with either population.
The pooled row shows $H_o = 0.375 \ll H_e = 0.494$, an apparent heterozygote deficit with $F \approx 0.24$ — out of nowhere, if you didn’t know the sample was structured.
The pooled $H_e$ (0.494) is also higher than either population’s own $H_e$ (0.320, 0.426). Pooling mixes allele frequencies toward the middle ($p = 0.553$, near 0.5 where $H_e$ peaks), so the expected diversity of the mixture exceeds the diversity of any real component.

This is the Wahlund effect: combining structured populations inflates apparent expected heterozygosity and produces an apparent heterozygote deficit — mimicking inbreeding even when every subpopulation is perfectly outbred.

Why does this matter for your work? If you compute heterozygosity from a mixed VCF without accounting for population structure, substructure masquerades as inbreeding. It also runs the other way: the PCA section below is exactly the tool you use to detect the structure that causes Wahlund effects — the connection between these two sections is the heart of “everything connects.”

In the red-spruce dataset, Population_Homozygosity computed across structured regions contains this same mixing effect, which is one reason it is analyzed by region and alongside PC1/PC2.

Exercises: heterozygosity

Exercise 6 — Ho and He for the two populations

Using pop_1 and pop_2 from Part II (N = 20 each):

Compute $H_o$ (observed frequency of Aa) and $H_e$ (from the observed $p$) for each population.
Compute $H_o$ and $H_e$ for the pooled population.
Is the pooled $H_o$ below the pooled $H_e$? Is the deficit larger or smaller than in the N = 500 demonstration above? Why?

Solution

Ho_pop1 <- sum(pop_1$genotype == "Aa") / nrow(pop_1)
Ho_pop2 <- sum(pop_2$genotype == "Aa") / nrow(pop_2)
Ho_pool <- sum(pooled_population$genotype == "Aa") / nrow(pooled_population)

He_pop1 <- 2 * p_A_1 * q_a_1
He_pop2 <- 2 * p_A_2 * q_a_2
He_pool <- 2 * p_A_p * q_a_p

data.frame(
  group = c("Pop_1", "Pop_2", "Pooled"),
  Ho = round(c(Ho_pop1, Ho_pop2, Ho_pool), 3),
  He = round(c(He_pop1, He_pop2, He_pool), 3),
  F = round(1 - c(Ho_pop1, Ho_pop2, Ho_pool) /
              c(He_pop1, He_pop2, He_pool), 3)
)

   group  Ho    He      F
1  Pop_1 0.5 0.420 -0.190
2  Pop_2 0.5 0.455 -0.099
3 Pooled 0.5 0.499 -0.003

With seed 42, all three groups happen to show $H_o \approx 0.50$ (and the within-pop $H_e$ values are 0.420 and 0.455), so the Wahlund deficit is essentially invisible: $F_{pool} \approx -0.003$ — a hair of heterozygote excess — even though the populations are just as differentiated as in the big simulation.

The reason is noise swamping signal: with N = 20, the standard error of a proportion near 0.5 is $\sqrt{0.25/20} \approx 0.11$, while the expected Wahlund deficit for this configuration is $H_{e,pool} - \overline{H_{o,within}} \approx 0.49 - 0.37 \approx 0.12$ — the same order as the noise. With N = 500, the standard error shrinks to $\approx 0.022$ and the deficit stands out clearly.

Practical takeaway: heterozygote deficits (and excesses) in small samples are uninterpretable without either larger N or replicated loci. This is why genome-scale analyses average over thousands of loci — the per-locus noise cancels while a systematic deficit (structure, inbreeding) persists.

Exercise 7 — Average across loci

For the geno matrix from Part I:

Add a per-locus inbreeding coefficient $F = 1 - H_o/H_e$ to het_table.
Which single locus departs most from H–W in each direction?
Compute the mean $H_o$, mean $H_e$, and multi-locus $F$ across all 10 loci. Why is averaging across loci more trustworthy than any single locus?

Solution

# Compute F from the unrounded Ho and He, then attach it to the table
het_table$F <- round(1 - Ho / He, 3)
het_table

            locus    Ho    He   diff      F
locus_1   locus_1 0.267 0.480 -0.213  0.444
locus_2   locus_2 0.500 0.473  0.027 -0.058
locus_3   locus_3 0.367 0.473 -0.106  0.224
locus_4   locus_4 0.333 0.464 -0.131  0.282
locus_5   locus_5 0.467 0.480 -0.013  0.028
locus_6   locus_6 0.433 0.433  0.001 -0.001
locus_7   locus_7 0.367 0.486 -0.119  0.246
locus_8   locus_8 0.400 0.480 -0.080  0.167
locus_9   locus_9 0.333 0.464 -0.131  0.282
locus_10 locus_10 0.567 0.495  0.072 -0.145

mean(het_table$Ho)

[1] 0.4034

mean(het_table$He)

[1] 0.4728

1 - mean(het_table$Ho) / mean(het_table$He)

[1] 0.1467851

With seed 123: the strongest deficit is locus_1 ($F = 0.44$), the strongest excess is locus_10 ($F = -0.15$), and the multi-locus value is $F \approx 0.147$.

Averaging across loci works because each locus is an independent draw of the same underlying process. Random per-locus deviations — which push individual $F$ values both up and down — cancel out, while a systematic process (inbreeding, structure, or in our case the simulation’s genotype probabilities) shifts all loci in the same direction and survives the average. This is the same logic as Exercise 6 writ large: more data, same signal, less noise. It is also exactly why real genomic studies report genome-wide average $H_o/H_e$ rather than single-locus values.

Exercise 8 — Connect to the real data

Using the red-spruce dataset:

Plot Genetic_Diversity against Population_Homozygosity.
Fit lm(Population_Homozygosity ~ Genetic_Diversity, data = traits).
Interpret the slope in light of the identity $H_e + \text{homozygosity} = 1$ from the worked example above.

Solution sketch

plot(traits$Genetic_Diversity, traits$Population_Homozygosity,
     pch = 19, xlab = "Genetic diversity", ylab = "Population homozygosity")
summary(lm(Population_Homozygosity ~ Genetic_Diversity, data = traits))

For a single biallelic locus, expected homozygosity is $p^2 + q^2 = 1 - 2pq = 1 - H_e$, so the two quantities are deterministically complementary. If these columns were computed from the same set of loci and the same weighting, you should find a strongly negative relationship — though not necessarily a slope of exactly $-1$, since the columns may use different locus sets, weightings, or estimators. Scatter around a strong negative trend is expected with real, noisy, structured data; deviations can also reflect population structure (Wahlund again) inflating homozygosity in mixed samples. Check the outlier populations against their Region on the PCA plot.

PCA for population structure

On the simulated matrix:

pca_sim <- prcomp(geno, scale. = TRUE)
plot(pca_sim$x[, 1], pca_sim$x[, 2],
     xlab = "PC1", ylab = "PC2", main = "Simulated genotype PCA")

Now look at the real PC1/PC2 already computed for the red-spruce trees (from actual genetic markers, not simulated):

traits <- read.table("data/raw/red_spruce_fitness_traits.txt", header = TRUE, sep = "\t")

traits$Region <- factor(traits$Region)
plot(traits$PC1, traits$PC2,
     col = as.integer(traits$Region), pch = 19,
     xlab = "PC1", ylab = "PC2", main = "Genetic PCA colored by region")
legend("topright", legend = levels(traits$Region),
       col = 1:nlevels(traits$Region), pch = 19)

In population genetics, PC1/PC2 from a genetic marker PCA usually capture the strongest axes of population structure — often correlated with geography (isolation by distance) or major environmental gradients.

Worked example: what a PCA of the two-population simulation would show

Go back to the Wahlund simulation (two populations, $p = 0.8$ vs. $p = 0.3$, N = 500 each) and imagine running prcomp() on a 1000 $\times$ 1 matrix of genotypes. Because there is only one locus, the “PCA” is trivial — but if you extend the simulation to hundreds of loci that all differ between the populations in the same direction (like locally adapted loci), PC1 would capture exactly one thing: which population an individual belongs to.

That is what PC1/PC2 do in real data. When the red-spruce PCA shows points colored by Region separating along PC1, you are looking at the multilocus version of the allele-frequency difference $p_{A,1} = 0.70$ vs. $p_{A,2} = 0.35$ from Part II. The PCA doesn’t create information — it summarizes thousands of per-locus frequency differences into two axes you can plot.

The loop closes here: structure (Part II pooling) → heterozygote deficit (Wahlund) → PCA detects that structure → and the GEA/RDA step must account for it, because a climate variable correlated with region would otherwise look like adaptation when it is really just ancestry. This is why PC1 sits in the red-spruce dataset alongside the climate variables.

How everything connects

Everything in this module is one pipeline applied at increasing scale. Trace it end-to-end with a single question: “why is this locus heterozygote-deficient?”

Step	Concept	Where you computed it	Red-spruce equivalent
1	Count genotypes	`table(pop_1$genotype)`	—
2	Allele frequencies $p, q$	`colMeans(geno)/2`, allele counting	—
3	H–W expectations $p^2, 2pq, q^2$	`prob = c(p^2, 2pq, q^2)`	—
4	Observed vs. expected heterozygosity	`Ho`, `He`, `het_table`	`Genetic_Diversity`
5	The gap: $F = 1 - H_o/H_e$	Exercise 7	`Family_Homozygosity` (conceptually related)
6	Structure causes the gap	Wahlund effect (N = 500 demo)	`Region`, `PC1`, `PC2`
7	Structure detected multivariately	`prcomp()` PCA plot	`PC1`, `PC2` columns
8	Structure + environment	`cor.test()`, `lm()`	`Elevation`, `Latitude`
9	Many loci + environment jointly	`vegan::rda()` (stretch)	RDA / genomic offset

Three threads run through the whole table:

The same two quantities everywhere. Every row is computed from genotype counts and sample size. Heterozygosity, F, the Wahlund effect, PCA loadings, and even the RDA ordination are all functions of allele frequencies — the only thing that changes is how many loci and how many environmental axes you let in.
Structure is the recurring confounder. Pooling populations inflates apparent diversity (Part II, Wahlund), heterozygote deficits can mean either inbreeding or structure (Part III), and the PCA that reveals structure (PC1/PC2) is also the covariate you must control before interpreting any genotype–environment association (Part III’s lm() and rda()). A climate correlation that survives after accounting for structure is evidence of local adaptation; one that disappears was ancestry all along.
Scale changes confidence, not logic. With N = 20 at one locus, noise dominates (Exercise 6). With N = 500, the Wahlund effect is obvious. With 10 loci, averaging helps (Exercise 7). With genome-wide SNPs, per-locus estimates are noise but multivariate summaries (PCA, RDA) become powerful. The math you did by hand on a 2-allele, 20-individual toy is the same math inside adegenet, LEA, and your genomic-offset pipelines.

When you reach the stretch exercise — rda() with Genetic_Diversity, Genetic_Load, PC1, PC2 as responses and Elevation, Latitude as predictors — notice what it is: rows 8 and 9 of the table. A correlation tests one response–predictor pair; the RDA tests all of them jointly, in the same way the PCA summarized all loci jointly. Multivariate = the same idea, more axes.

Roadmap beyond this module

When you’re ready to work with real genotype data instead of pre-computed summaries:

vcfR::read.vcfR() reads a VCF file into R
adegenet::df2genind() / vcfR::vcfR2genind() convert genotypes into population-genetics-aware objects with built-in allele-frequency, heterozygosity, and Fst functions
LEA (Bioconductor) and vegan::rda() scale the GEA approach above to genome-wide SNP data

None of this is required for this module’s exercises — it’s here so you know where basics you learn connects to your actual genomic analysis pipelines.

Wrap Up

Idea check

Answer these before the practical exercises.

What does an alternate allele frequency of 0.3 at a locus mean biologically?
Why can observed heterozygosity (Ho) be lower than expected heterozygosity (He) in a real population?
What do you think Population_Homozygosity is measuring, and how would you expect it to relate to Genetic_Diversity — positively or negatively?
In a PCA built from genetic markers, what do PC1 and PC2 typically represent in a population-genetics context?
Why test cor.test(Genetic_Diversity, Elevation) rather than just assuming a relationship exists because both vary by population?
What is the conceptual difference between running a separate lm() for each genetic variable versus running one rda() with all of them as a response matrix?

Practical exercises

Build the simulated genotype matrix (30 individuals × 10 loci) exactly as shown, and compute the alternate allele frequency for every locus.
Compute observed and expected heterozygosity for each locus and identify which locus shows the largest Ho − He gap.
Run prcomp() on the simulated genotype matrix and plot PC1 vs. PC2.
Load the real red-spruce dataset and plot the real PC1 vs. PC2, colored by Region. Does the pattern look like discrete clusters or a continuous gradient?
Run cor.test() between Population_Homozygosity and Elevation, and separately between Genetic_Diversity and Latitude. Report the correlation coefficient, p-value, and a one-sentence interpretation for each.
Fit lm(Genetic_Load ~ Elevation + Latitude, data = traits) and interpret the two coefficients — which predictor has a stronger association, and in which direction?
Stretch exercise: run the vegan::rda() example above using Genetic_Diversity, Genetic_Load, PC1, PC2 as responses and Elevation, Latitude as predictors. Produce the ordination biplot and describe what you see in 2–3 sentences.
Save every summary table (allele frequencies, heterozygosity table, correlation results) and every plot into results/.

Solutions for the simulation-based items (1–3) are worked through in Part I, the Part II exercises, and Exercise 7 — attempt them before unfolding those callouts.

Homework

Build src/02_pop_genetics_intro.R, combining the simulated-data fundamentals with a real-data GEA-style analysis.

Requirements:

Part A (fundamentals): allele frequency and Ho/He table for the simulated genotype matrix, plus the PCA plot.
Part B (applied): using the real dataset, test associations between at least two genetic-parameter columns (e.g., Population_Homozygosity, Genetic_Diversity, Genetic_Load) and at least two environmental predictors (Elevation, Latitude, Longitude), using correlation and/or lm().
A written interpretation (half a page, as comments or a results/README.md) answering: Do the genetic parameters show a spatial pattern with latitude or elevation? What ecological or demographic process might explain this (e.g., isolation by distance, local adaptation, drift in small/marginal populations)?
Bonus (optional): complete the vegan::rda() ordination and include the biplot with a short interpretation of which environmental axis drives more separation among genetic variables.

Submission format: src/02_pop_genetics_intro.R, all generated outputs in results/, and the written interpretation.