tests/testthat/test-geneflow.R

In bodkan/slendr: A Simulation Framework for Spatiotemporal Population Genetics

test_that("non-overlapping geneflow leads to error", {
  map <- readRDS("map.rds")
  pop1 <- population("pop1", N = 100, time = 10000,
                      center = c(0, 10), radius = 1, map = map, intersect = FALSE)
  pop2 <- population("pop2", N = 100, time = 10000,
                     center = c(0, -10), radius = 1, map = map, intersect = FALSE)
  expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = 0.1),
               "No overlap between population ranges of pop1 and pop2 at time 1000.")
})

test_that("non-overlapping geneflow passes if the check is explicitly turned off", {
  map <- readRDS("map.rds")
  pop1 <- population("pop1", N = 100, time = 10000,
                     center = c(0, 10), radius = 1, map = map, intersect = FALSE)
  pop2 <- population("pop2", N = 100, time = 10000,
                     center = c(0, -10), radius = 1, map = map, intersect = FALSE)
  expect_s3_class(
    gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = 0.1, overlap = FALSE),
    "data.frame"
  )
})

test_that("populations must be present within the given gene-flow time window (spatial)", {
  map <- readRDS("map.rds")
  pop1 <- population("pop1", N = 100, time = 100,
                      center = c(0, 10), radius = 1, map = map, intersect = FALSE)
  pop2 <- population("pop2", N = 100, time = 100,
                     center = c(0, -10), radius = 1, map = map, intersect = FALSE)
  expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = 0.1),
               "Both .* and .* must be already present within the gene-flow window \\d+-\\d+")
})

test_that("populations must be present within the given gene-flow time window (non-spatial)", {
  pop1 <- population("pop1", N = 100, time = 100)
  pop2 <- population("pop2", N = 100, time = 100)
  expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = 0.1),
               "Both .* and .* must be already present within the gene-flow window \\d+-\\d+")
})

test_that("gene-flow rate must be a value between 0 and 1", {
  pop1 <- population("pop1", N = 100, time = 100)
  pop2 <- population("pop2", N = 100, time = 100)

  error_msg <- "Gene-flow rate must be a numeric value between 0 and 1"
  expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = -0.1), error_msg)
  expect_error(gene_flow(from = pop1, to = pop2, start = 1000, end = 0, rate = 25), error_msg)
})

test_that("populations must be already created for a gene flow to happen (forward model)", {
  skip_if(!check_dependencies(python = TRUE))
  init_env(quiet = TRUE)

  pop1 <- population("pop1", N = 100, time = 100)
  pop2 <- population("pop2", N = 100, time = 200, parent = pop1)

  error_msg <- "^Both .* and .* must be already present within the gene-flow window .*-.*$"
  # neither population missing:
  expect_error(gene_flow(from = pop1, to = pop2, start = 50, end = 80, rate = 0.1), error_msg)
  # pop1 present, pop2 not created yet:
  expect_error(gene_flow(from = pop1, to = pop2, start = 50, end = 130, rate = 0.1), error_msg)
  expect_error(gene_flow(from = pop1, to = pop2, start = 200, end = 230, rate = 0.1), error_msg)
  # (both present if gene flow starts one generation later)
  expect_silent(gene_flow(from = pop1, to = pop2, start = 201, end = 230, rate = 0.1))
  # pop1 not yet created
  expect_error(gene_flow(from = pop1, to = pop2, start = 100, end = 230, rate = 0.1), error_msg)

  # full run of a scenario which should pass gene_flow(), compile_model() and create a tree sequence
  expect_silent({
    gf <- gene_flow(from = pop1, to = pop2, start = 201, end = 230, rate = 0.1)
    model <- compile_model(list(pop1, pop2), generation_time = 1, simulation_length = 201)
  })
  expect_s3_class(ts_msprime <- msprime(model, sequence_length = 100, recombination_rate = 0),
                  "slendr_ts")
  expect_s3_class(ts_slim <- slim(model, sequence_length = 100, recombination_rate = 0),
                  "slendr_ts")
})

# possibly reduntant given the above, but a useful test of an issue reported on Github
# (for quite a long time, gene_flow() failed correctly but with confusing error message
# due to weird code duplication inside gene_flow()) -- this test is added to check
# that the issue has been solved
# https://github.com/bodkan/slendr/issues/132
test_that("gene_flow() behaves as expected in some concrete situations", {
  anc <- population("ancestor", time = 9e6, remove=6e6, N=100)
  popA <- population("A", time = 8e6, parent=anc, N=100)
  popB <- population("B", time = 7e6, parent=anc, N=100)

  # this used to say: "Both A and B must be present within the gene-flow window 4e+06-4000005"
  expect_error(
    gene_flow(from = popA, to = popB, start = 4e6, end = (4e6+5), rate = 0.3),
    "Inconsistent time direction implied by populations and the gene flow event"
  )

  # this used to say: "Specified times are not consistent with the assumed direction of time
  # (gene flow A -> B in the time window 8e+06-7999995)"
  expect_error(
    gene_flow(from = popA, to = popB, start = 8e6, end = (8e6-5), rate = 0.3),
    "Both A and B must be already present within the gene-flow window .*-.*"
  )
})

test_that("gene flow must take longer than a generation (forward models)", {
  init_env(quiet = TRUE)

  p1 <- population("p1", time = 1, N = 100)
  p2 <- population("p2", time = 100, N = 100, parent = p1)
  p3 <- population("p3", time = 200, N = 100, parent = p2)

  # gene flow before the existence of either population
  expect_error(
    gf <- gene_flow(from = p2, to = p3, rate = 0.5, start = 10, end = 80),
    "Both p2 and p3 must be already present within the gene-flow window 10-80"
  )
  # gene flow before the existence of the target population (off by 1 "year")
  expect_error(
    gf <- gene_flow(from = p2, to = p3, rate = 0.5, start = 200, end = 250),
    "Both p2 and p3 must be already present within the gene-flow window 200-250"
  )
  # shifting the start by 1 "year" fixes the error
  expect_s3_class(
    gf <- gene_flow(from = p2, to = p3, rate = 0.5, start = 201, end = 250),
    "data.frame"
  )

  # problem: gene flow can be shorter than the span of a generation

  # single invalid gene flow
  gf <- gene_flow(from = p3, to = p1, rate = 0.5, start = 221, end = 230)
  expect_error(
    model <- compile_model(list(p1, p2, p3), generation_time = 20, gene_flow = gf, simulation_length = 1000),
    "Gene flows cannot be shorter than the generation time of 20"
  )

  # multiple invalid gene flows
  gf <- list(
    gene_flow(from = p2, to = p3, rate = 0.5, start = 201, end = 210),
    gene_flow(from = p3, to = p1, rate = 0.5, start = 221, end = 230)
  )

  expect_error(
    model <- compile_model(list(p1, p2, p3), generation_time = 20,
                           gene_flow = gf, simulation_length = 1000),
    "Gene flows cannot be shorter than the generation time of 20"
  )

  # single invalid gene flows among multiple
  gf <- list(
    gene_flow(from = p2, to = p3, rate = 0.5, start = 201, end = 250),
    gene_flow(from = p3, to = p1, rate = 0.5, start = 221, end = 230)
  )

  expect_error(
    model <- compile_model(list(p1, p2, p3), generation_time = 20,
                           gene_flow = gf, simulation_length = 1000),
    "Gene flows cannot be shorter than the generation time of 20"
  )

  # properly timed gene flows work OK
  t <- 20
  gf <- list(
    gene_flow(from = p2, to = p3, rate = 0.5, start = 300, end = 300 + t),
    gene_flow(from = p3, to = p1, rate = 0.5, start = 500, end = 500 + t)
  )
  expect_s3_class(
    model <- compile_model(list(p1, p2, p3), generation_time = t,
                           gene_flow = gf, simulation_length = 1000),
    "slendr_model"
  )
})

test_that("gene flow must take longer than a generation (backward models)", {
  init_env(quiet = TRUE)

  p1 <- population("p1", time = 1000, N = 100)
  p2 <- population("p2", time = 800, N = 100, parent = p1)
  p3 <- population("p3", time = 600, N = 100, parent = p2)

  # gene flow before the existence of either population
  expect_error(
    gf <- gene_flow(from = p2, to = p3, rate = 0.5, start = 1000, end = 900),
    "Both p2 and p3 must be already present within the gene-flow window 1000-900",
  )
  # gene flow before the existence of the target population (off by 1 "year")
  expect_error(
    gf <- gene_flow(from = p2, to = p3, rate = 0.5, start = 600, end = 500),
    "Both p2 and p3 must be already present within the gene-flow window 600-500"
  )
  # shifting the start by 1 "year" fixes the error
  expect_s3_class(
    gf <- gene_flow(from = p2, to = p3, rate = 0.5, start = 599, end = 500),
    "data.frame"
  )

  # problem: gene flow can be shorter than the span of a generation

  # single invalid gene flow
  gf <- gene_flow(from = p3, to = p1, rate = 0.5, start = 500, end = 490)
  expect_error(
    model <- compile_model(list(p1, p2, p3), generation_time = 20, gene_flow = gf, simulation_length = 1000),
    "Gene flows cannot be shorter than the generation time of 20"
  )

  # multiple invalid gene flows
  gf <- list(
    gene_flow(from = p2, to = p3, rate = 0.5, start = 500, end = 490),
    gene_flow(from = p3, to = p1, rate = 0.5, start = 100, end = 90)
  )
  expect_error(
    model <- compile_model(list(p1, p2, p3), generation_time = 20,
                           gene_flow = gf, simulation_length = 1000),
    "Gene flows cannot be shorter than the generation time of 20"
  )

  # single invalid gene flows among multiple
  gf <- list(
    gene_flow(from = p2, to = p3, rate = 0.5, start = 500, end = 490),
    gene_flow(from = p3, to = p1, rate = 0.5, start = 100, end = 90)
  )
  expect_error(
    model <- compile_model(list(p1, p2, p3), generation_time = 20,
                           gene_flow = gf, simulation_length = 1000),
    "Gene flows cannot be shorter than the generation time of 20"
  )

  # properly timed gene flows work OK
  t <- 20
  gf <- list(
    gene_flow(from = p2, to = p3, rate = 0.5, start = 500, end = 500 - t),
    gene_flow(from = p3, to = p1, rate = 0.5, start = 100, end = 100 - t)
  )
  expect_s3_class(
    model <- compile_model(list(p1, p2, p3), generation_time = t,
                           gene_flow = gf, simulation_length = 1000),
    "slendr_model"
  )
})