5  Genetic drift

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.3.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(janitor)

Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

N<- 10
population <- tibble(ID = 1:(2*N), allele = c(rep("blue", N), rep("red",N)), value=c(rep(0,N),rep(1,N)))

generation<-0
new_population<-slice_sample(population,n=nrow(population), replace = TRUE)
n_red <- sum(new_population$value)
while (n_red>0 & n_red<2*N){
  generation <- generation + 1
  n_red <- sum(new_population$value)
  print(paste(generation, n_red))
  new_population<-slice_sample(new_population,n=nrow(population), replace = TRUE)
}

[1] "1 9"
[1] "2 10"
[1] "3 13"
[1] "4 12"
[1] "5 13"
[1] "6 13"
[1] "7 13"
[1] "8 17"
[1] "9 15"
[1] "10 14"
[1] "11 16"
[1] "12 17"
[1] "13 15"
[1] "14 16"
[1] "15 18"
[1] "16 14"
[1] "17 16"
[1] "18 14"
[1] "19 12"
[1] "20 13"
[1] "21 13"
[1] "22 12"
[1] "23 10"
[1] "24 14"
[1] "25 16"
[1] "26 18"
[1] "27 17"
[1] "28 18"
[1] "29 18"
[1] "30 18"
[1] "31 16"
[1] "32 18"
[1] "33 17"
[1] "34 19"
[1] "35 19"
[1] "36 19"
[1] "37 20"

5.1 Binomial version

# Initial state
npops <- 5
pop_size <- 100
total_alleles <- 2 * pop_size
p <- 0.5 # Frequency of allele A
generations <- 300

ps <- matrix(rep(0, npops*generations), ncol = npops)
i <-0
for (pop in 1:npops){
  p <-0.5
  for (i in 1:generations){
      # Randomly draw 2N alleles based on frequency p
      # Each draw is binomial sampling
      num_A <- rbinom(1,total_alleles, p)
      
      # Update frequency
      p <- num_A / total_alleles
      ps[i,pop] <- p
      
      # print(f"Gen {gen+1}: Freq A = {p}")
      
      # If allele is fixed or lost, stop
  }
}

ps_df <- as.tibble(ps)

Warning: `as.tibble()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` instead.
ℹ The signature and semantics have changed, see `?as_tibble`.

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.
ℹ The deprecated feature was likely used in the tibble package.
  Please report the issue at <https://github.com/tidyverse/tibble/issues>.

names(ps_df) <- LETTERS[1:npops] 
ps_df <-   mutate(ps_df, generation = 1:generations) |>
  select(generation,names(ps_df))

ps_df |>
  pivot_longer(A:E,names_to= "population", values_to = "allele_frequency") |>
  ggplot(aes(x = generation, y = allele_frequency, colour = population)) +
    geom_point(size = 1) +
    geom_line(linewidth = 0.3) +
    labs(x = "Generation",
         y = "Allele Frequency") +
    theme_minimal()