library(gstudio)
library(dplyr)11 Simulation & Parentage
This chapter covers the forward-time population simulation framework and parentage analysis tools in gstudio.
11.1 Building Populations from Frequencies
11.1.1 Allele Frequency Specification
The simulation tools start with specifying allele frequencies per population and locus. The make_populations() function requires a data.frame with Population, Locus, Allele, and Frequency columns:
freqs <- expand.grid(Allele = c("01", "02"),
Locus = c("Loc1", "Loc2", "Loc3"),
Population = c("A", "B", "C"),
stringsAsFactors = FALSE) |>
mutate(Frequency = c(0.2, 0.8, 0.3, 0.7, 0.4, 0.6, # A
0.7, 0.3, 0.5, 0.5, 0.6, 0.4, # B
0.4, 0.6, 0.3, 0.7, 0.2, 0.8), # C
Population = factor(Population, ordered = TRUE)) |>
select(Population, Locus, Allele, Frequency)11.1.2 Creating Populations
pop <- freqs |>
make_populations(N = 30)
pop |>
head() Population ID Loc1 Loc2 Loc3
1 A 1 02:02 01:01 01:02
2 A 2 01:02 02:02 02:02
3 A 3 01:02 01:02 02:02
4 A 4 01:02 02:02 02:02
5 A 5 02:02 02:02 01:02
6 A 6 01:01 02:02 02:02pop |>
count(Population) Population n
1 A 30
2 B 30
3 C 30You can specify an inbreeding level with F:
pop_inbred <- freqs |>
make_populations(N = 30, F = 0.3)11.2 Migration Models
The migration_matrix() function creates a K x K migration matrix where rows sum to 1. Element [i,j] is the proportion of population i that originates from population j.
11.2.1 Island Model
Equal migration among all populations:
migration_matrix(3, model = "island", m = 0.05) Pop1 Pop2 Pop3
Pop1 0.950 0.025 0.025
Pop2 0.025 0.950 0.025
Pop3 0.025 0.025 0.950A count names the populations Pop1, Pop2, …. For simulations, pass the population names instead (as below), so the matrix matches the Population column that simulate_pop() migrates between.
11.2.2 1D Stepping Stone
Migration only between adjacent populations:
c("A", "B", "C", "D") |>
migration_matrix(model = "stepping_stone_1d", m = 0.1) A B C D
A 0.90 0.10 0.00 0.00
B 0.05 0.90 0.05 0.00
C 0.00 0.05 0.90 0.05
D 0.00 0.00 0.10 0.9011.2.3 2D Stepping Stone
Grid-based adjacency:
migration_matrix(4, model = "stepping_stone_2d",
m = 0.1, nr = 2, nc = 2) Pop1 Pop2 Pop3 Pop4
Pop1 0.90 0.05 0.05 0.00
Pop2 0.05 0.90 0.00 0.05
Pop3 0.05 0.00 0.90 0.05
Pop4 0.00 0.05 0.05 0.9011.2.4 Distance-Based
Migration rate decays with distance:
coords <- matrix(c(0, 0, 1, 0, 0.5, 1), ncol = 2, byrow = TRUE)
migration_matrix(3, model = "distance", m = 0.1, coords = coords) Pop1 Pop2 Pop3
Pop1 0.9000000 0.0527864 0.0472136
Pop2 0.0527864 0.9000000 0.0472136
Pop3 0.0500000 0.0500000 0.900000011.3 Migration Events
A migration_event pairs a migration matrix with a generation interval, allowing temporal regime changes:
# Low migration for generations 1-50
evt1 <- c("A", "B", "C") |>
migration_matrix(model = "island", m = 0.05) |>
migration_event(start = 1, end = 50)
# High migration from generation 51 onward
evt2 <- c("A", "B", "C") |>
migration_matrix(model = "island", m = 0.20) |>
migration_event(start = 51)11.4 Mutation Models
The mutation_model() function defines how alleles mutate during simulation:
# Infinite Alleles Model
mm_iam <- mutation_model(rate = 0.001, model = "iam")
# K-Allele Model (10 possible states)
mm_kam <- mutation_model(rate = 0.001, model = "kam", k = 10)
# Stepwise Mutation Model
mm_smm <- mutation_model(rate = 0.001, model = "smm")11.5 Forward-Time Simulation
The simulate_pop() function runs a multi-generation forward-time simulation. Each generation cycle: census, migrate, mate, mutate.
11.5.1 Basic Simulation
result <- pop |>
simulate_pop(ngen = 20, verbose = FALSE)
result |>
head() ID Population Loc1 Loc2 Loc3
1 1 A 01:02 01:01 01:02
2 2 A 02:02 02:02 02:02
3 3 A 02:02 02:02 02:02
4 4 A 02:02 02:02 01:02
5 5 A 01:02 02:02 02:02
6 6 A 02:02 01:02 02:0211.5.2 With Migration
evt <- pop |>
pull(Population) |>
unique() |>
migration_matrix(model = "island", m = 0.05) |>
migration_event(start = 1)
result <- pop |>
simulate_pop(ngen = 20, migration = evt, verbose = FALSE)11.5.3 With Mutation
mm <- mutation_model(rate = 0.01, model = "iam")
result <- pop |>
simulate_pop(ngen = 20, mutation = mm, verbose = FALSE)11.5.4 With Mixed Mating (Selfing)
result <- pop |>
simulate_pop(ngen = 20, selfing_rate = 0.3, verbose = FALSE)11.5.5 Census Snapshots
Save population snapshots at regular intervals:
result <- pop |>
simulate_pop(ngen = 100, census_interval = 25, census_dir = tempdir(), verbose = FALSE)
# Snapshots saved as gen_0025.rds, gen_0050.rds, etc.11.5.6 Complete Worked Example
Simulate drift and migration in a 3-population system, then analyze the outcome:
# Create starting populations
pop <- freqs |>
make_populations(N = 50)
# Set up island model migration among the simulated populations
evt <- pop |>
pull(Population) |>
unique() |>
migration_matrix(model = "island", m = 0.02) |>
migration_event(start = 1)
# Run simulation
set.seed(42)
final <- pop |>
simulate_pop(ngen = 50, migration = evt, verbose = FALSE)
# Analyze structure in the final generation
final |>
genetic_structure(stratum = "Population", mode = "Gst") Locus Gst Hs Ht p.value
1 Loc1 0.03860354 0.4822222 0.5015852 NA
2 Loc2 0.02426476 0.4580471 0.4694379 NA
3 Loc3 0.11745245 0.3934007 0.4457558 NA
4 Multilocus 0.05866043 1.3336700 1.4167789 NA# Compare diversity before and after
list(Start = pop, Final = final) |>
lapply(genetic_diversity, stratum = "Population", mode = "He") |>
bind_rows(.id = "Generation") |>
group_by(Generation, Stratum) |>
summarize(mean_He = mean(He), .groups = "drop")# A tibble: 6 × 3
Generation Stratum mean_He
<chr> <chr> <dbl>
1 Final A 0.485
2 Final B 0.470
3 Final C 0.365
4 Start A 0.416
5 Start B 0.460
6 Start C 0.39111.6 Parentage Analysis
11.6.1 Mating and Offspring
The mate() function creates offspring from two parents:
f <- data.frame(
Locus = rep(paste0("Loc", 1:4), each = 3),
Allele = rep(LETTERS[1:3], times = 4),
Frequency = rep(c(1/3, 1/3, 1/3), times = 4)
)
adults <- f |>
make_population(N = 20)
offs <- adults |>
slice(1) |>
mate(slice(adults, 2), N = 5)
offs ID Loc1 Loc2 Loc3 Loc4
1 1 A:C A:C A:B B:C
2 1 A:C A:C A:A B:C
3 1 A:C A:C A:B B:C
4 1 A:C A:C A:A A:B
5 1 B:C A:C A:B A:B11.6.2 Maternal Subtraction
The minus_mom() function removes the maternal contribution from offspring genotypes. It takes a single data.frame where mothers have OffID = 0 and offspring have OffID != 0:
# Combine mother (OffID = 0) and offspring into one data.frame
family <- bind_rows(adults |>
slice(1) |>
mutate(OffID = 0),
offs |>
mutate(OffID = row_number()))
family |>
minus_mom() ID Loc1 Loc2 Loc3 Loc4 OffID
2 1 C A A:B B 1
3 1 C A A B 2
4 1 C A A:B B 3
5 1 C A A B 4
6 1 C A A:B B 511.6.3 Fractional Paternity
The paternity() function estimates the probability of each potential father siring each offspring:
offs <- offs |>
mutate(OffID = row_number(),
MomID = adults$ID[1])
pat <- offs |>
paternity(slice(adults, 1), adults)
pat |>
head() MomID OffID DadID Fij
1 1 1 2 0.5714286
2 1 1 12 0.1428571
3 1 1 16 0.1428571
4 1 1 17 0.1428571
5 1 2 2 0.6153846
6 1 2 16 0.1538462The Fij column gives the fractional paternity probability. Use strict = TRUE to retain only unambiguous assignments:
offs |>
paternity(slice(adults, 1), adults, strict = TRUE)[1] MomID OffID DadID Fij
<0 rows> (or 0-length row.names)11.6.4 Transition Probabilities and the Paternity Index
transition_probability() gives the Mendelian probability of an offspring genotype from a pair of parents, as a product across loci (multilocus = FALSE returns the per-locus values):
mom <- slice(adults, 1)
dad <- slice(adults, 2)
transition_probability(slice(offs, 1), mom, dad)[1] 0.0625transition_probability(slice(offs, 1), mom, dad, multilocus = FALSE)Loc1 Loc2 Loc3 Loc4
0.5 0.5 0.5 0.5 To judge a putative father, compare that value with the probability of the offspring given the mother and a random male. Pass only the mother and an allele frequency table from frequencies(); at each locus this is the chance the mother transmits an allele the offspring carries times the frequency of the offspring’s other allele. The ratio is the paternity index, how many times more likely the male is to be the father than a random male:
freqs <- adults |>
frequencies()
T_random <- transition_probability(slice(offs, 1), mom, freqs = freqs)
transition_probability(slice(offs, 1), mom, dad) / T_random[1] 22.7047Without freqs, the one-parent value is only a compatibility score (zero when the parent is excluded), not a likelihood, so do not use it as the denominator of a paternity index.
11.6.5 Single-Parent Exclusion
The parent_finder() function identifies compatible parents for each offspring using exclusion. It requires a data.frame with ID and OffID columns, where adults have OffID = 0:
# Prepare a combined data.frame: adults (OffID = 0) then offspring
locus_cols <- adults |>
column_class("locus")
combined <- bind_rows(adults |>
select(ID, all_of(locus_cols)) |>
mutate(OffID = 0),
offs |>
select(ID, all_of(locus_cols)) |>
mutate(OffID = row_number()))
combined |>
parent_finder() |>
head() ID OffID ParentID T
5 1 1 12 0.062500
2 1 1 2 0.031250
8 1 1 16 0.031250
9 1 1 17 0.031250
1 1 1 1 0.015625
6 1 1 13 0.015625Its T column is the one-parent compatibility score from transition_probability() (without freqs): it ranks the compatible parents of each offspring, but it is not a likelihood.
11.6.6 Identifying Bad Parents
The bad_parents() function flags parent-offspring pairs that are genetically incompatible:
offs |>
bad_parents(adults)11.7 Bundled Parentage Datasets
The cornus and cornus_florida datasets contain Cornus florida genotype data for parentage analysis:
data(cornus)
cornus |>
head() Population SampleID X.Coordinate Y.Coordinate Cf.G8 Cf.H18 Cf.N5 Cf.N10
1 2 203 3040 6146 157:165 117:119 170:170 189:191
2 2 204 3070 6114 159:163 101:103 166:170 193:201
3 2 205 3096 6148 151:157 113:119 193:195
4 2 206 3124 6136 153:161 117:119 197:201
5 2 207 3082 6126 159:163 105:119 158:170 193:193
6 2 208 3148 5180 181:181 119:119 170:170 199:201
Cf.O5
1 198:198
2 178:182
3 178:182
4 178:198
5 178:196
6 176:182data(cornus_florida)
cornus_florida |>
head() ID OffID X Y G8 H18 N5 N10 O5
1 226 0 1392 3534 162:180 114:114 124:126 192:192 185:195
2 232 0 1656 3414 158:180 112:112 124:126 184:192 185:185
3 234 0 1718 3330 158:180 112:96 128:128 184:192 185:185
4 300 0 1175 3114 180:188 112:116 126:126 198:200 191:195
5 305 0 1529 3237 154:170 122:124 124:126 188:192 181:195
6 432 0 1336 2748 164:180 114:116 124:126 198:202 185:19311.8 Without Pipes
pop <- make_populations(freqs, N = 50)
m <- migration_matrix(unique(pop$Population), model = "island", m = 0.02)
evt <- migration_event(m, start = 1)
final <- simulate_pop(pop, ngen = 50, migration = evt, verbose = FALSE)
genetic_structure(final, stratum = "Population", mode = "Gst")
adults <- make_population(f, N = 20)
offs <- mate(adults[1, ], adults[2, ], N = 5)
pat <- paternity(offs, adults[1, ], adults)