Package 'slendr'

Description

Animate the simulated population dynamics

Usage

animate_model(model, file, steps, gif = NULL, width = 800, height = 560)

Arguments

Value

If gif = NULL, return gganimate animation object. Otherwise a GIF file is saved and no value is returned.

Description

Calculate the area covered by the given slendr object

Usage

area(x)

Arguments

Value

Area covered by the input object. If a slendr_pop was given, a table with an population range area in each time point will be returned. If a slendr_region or slendr_world object was specified, the total area covered by this object's spatial boundary will be returned.

Examples

region_a <- region("A", center = c(20, 50), radius = 20)
region_b <- region("B", polygon = list(c(50, 40), c(70, 40), c(70, 60), c(50, 60)))
plot_map(region_a, region_b)

# note that area won't be *exactly* equal to pi*r^2:
#   https://stackoverflow.com/a/65280376
area(region_a)

area(region_b)

Description

Check that the required dependencies are available for slendr to work

Usage

check_dependencies(python = FALSE, slim = FALSE, quit = FALSE)

Arguments

Value

If quit = TRUE, no values is returned, if quit = FALSE, a scalar logical value is returned indicating whether or not the dependencies are present.

Description

This function inspects the Python environment which has been activated by the reticulate package and prints the versions of all slendr Python dependencies to the console.

Usage

check_env(verbose = TRUE)

Arguments

Value

Either TRUE (slendr Python environment is present) or FALSE (slendr Python environment is not present).

Examples

init_env()
check_env()

Description

Remove the automatically created slendr Python environment

Usage

clear_env(force = FALSE, all = FALSE)

Arguments

Value

No return value, called for side effects

Description

First, compiles the vectorized population spatial maps into a series of binary raster PNG files, which is the format that SLiM understands and uses it to define population boundaries. Then extracts the demographic model defined by the user (i.e. population divergences and gene flow events) into a series of tables which are later used by the built-in SLiM script to program the timing of simulation events.

Usage

compile_model(
  populations,
  generation_time,
  gene_flow = list(),
  direction = NULL,
  simulation_length = NULL,
  serialize = TRUE,
  path = NULL,
  overwrite = FALSE,
  force = FALSE,
  description = "",
  time_units = NULL,
  resolution = NULL,
  competition = NULL,
  mating = NULL,
  dispersal = NULL,
  extension = NULL
)

Arguments

Value

Compiled slendr_model model object which encapsulates all information about the specified model (which populations are involved, when and how much gene flow should occur, what is the spatial resolution of a map, and what spatial dispersal and mating parameters should be used in a SLiM simulation, if applicable)

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

Calculate the distance between a pair of spatial boundaries

Usage

distance(x, y, measure, time = NULL)

Arguments

Value

If the coordinate reference system was specified, a distance in projected units (i.e. meters) is returned. Otherwise the function returns a normal Euclidean distance.

Examples

# create two regions on a blank abstract landscape
region_a <- region("A", center = c(20, 50), radius = 20)
region_b <- region("B", center = c(80, 50), radius = 20)
plot_map(region_a, region_b)

# compute the distance between the centers of both population ranges
distance(region_a, region_b, measure = "center")

# compute the distance between the borders of both population ranges
distance(region_a, region_b, measure = "border")

Description

Expands the spatial population range by a specified distance in a given time-window

Usage

expand_range(
  pop,
  by,
  end,
  start,
  overlap = 0.8,
  snapshots = NULL,
  polygon = NULL,
  lock = FALSE,
  verbose = TRUE
)

Arguments

Details

Note that because slendr models have to accomodate both SLiM and msprime back ends, population sizes and times of events are rounded to the nearest integer value.

Value

Object of the class slendr_pop, which contains population parameters such as name, time of appearance in the simulation, parent population (if any), and its spatial parameters such as map and spatial boundary.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

Open an interactive browser of the spatial model

Usage

explore_model(model)

Arguments

Value

No return value, called in order to start an interactive browser-based interface to explore the dynamics of a slendr model

Description

This function extract a slendr model parameters used to compile a given model object or simulate a tree sequence

Usage

extract_parameters(data)

Arguments

Value

A list of data frames containing parameters of the model used when compiling a model object

Examples

init_env()

# load an example model and simulate a tree sequence from it
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))
ts <- msprime(model, sequence_length = 1e5, recombination_rate = 0)

# extract model parameters from a compiled model object as a list of data frames
extract_parameters(model)

# the function can also extract parameters of a model which simulated a
# tree sequence
extract_parameters(ts)

Description

Define a gene-flow event between two populations

Usage

gene_flow(from, to, rate, start, end, overlap = TRUE)

Arguments

Value

Object of the class data.frame containing parameters of the specified gene-flow event.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

Get the name of the current slendr Python environment

Usage

get_env()

Value

Name of the slendr Python environment

Description

This function attempts to activate a dedicated slendr Miniconda Python environment previously set up via setup_env.

Usage

init_env(quiet = FALSE)

Arguments

Value

No return value, called for side effects

Description

Merge two spatial slendr objects into one

Usage

join(x, y, name = NULL)

Arguments

Value

Object of the class slendr_region which encodes a standard spatial object of the class sf with several additional attributes (most importantly a corresponding slendr_map object, if applicable).

Examples

# create a blank abstract world 1000x1000 distance units in size
blank_map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# it is possible to construct custom landscapes (islands, corridors, etc.)
island1 <- region("island1", polygon = list(c(10, 30), c(50, 30), c(40, 50), c(0, 40)))
island2 <- region("island2", polygon = list(c(60, 60), c(80, 40), c(100, 60), c(80, 80)))
island3 <- region("island3", center = c(20, 80), radius = 10)
archipelago <- island1 %>% join(island2) %>% join(island3)

custom_map <- world(xrange = c(1, 100), c(1, 100), landscape = archipelago)

# real Earth landscapes can be defined using freely-available Natural Earth
# project data and with the possibility to specify an appropriate Coordinate
# Reference System, such as this example of a map of Europe

real_map <- world(xrange = c(-15, 40), yrange = c(30, 60), crs = "EPSG:3035")

Description

This function defines a displacement of a population along a given trajectory in a given time frame

Usage

move(
  pop,
  trajectory,
  end,
  start,
  overlap = 0.8,
  snapshots = NULL,
  verbose = TRUE
)

Arguments

Value

Object of the class slendr_pop, which contains population parameters such as name, time of appearance in the simulation, parent population (if any), and its spatial parameters such as map and spatial boundary.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

This function will execute a built-in msprime script and run a compiled slendr demographic model.

Usage

msprime(
  model,
  sequence_length,
  recombination_rate,
  samples = NULL,
  random_seed = NULL,
  verbose = FALSE,
  debug = FALSE,
  run = TRUE,
  path = NULL
)

Arguments

Value

A tree-sequence object loaded via Python-R reticulate interface function ts_read (internally represented by the Python object tskit.trees.TreeSequence). If the path argument was set, it will return the path as a single-element character vector.

Examples

init_env()

# load an example model
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# afr and eur objects would normally be created before slendr model compilation,
# but here we take them out of the model object already compiled for this
# example (in a standard slendr simulation pipeline, this wouldn't be necessary)
afr <- model$populations[["AFR"]]
eur <- model$populations[["EUR"]]
chimp <- model$populations[["CH"]]

# schedule the sampling of a couple of ancient and present-day individuals
# given model at 20 ky, 10 ky, 5ky ago and at present-day (time 0)
modern_samples <- schedule_sampling(model, times = 0, list(afr, 10), list(eur, 100), list(chimp, 1))
ancient_samples <- schedule_sampling(model, times = c(40000, 30000, 20000, 10000), list(eur, 1))

# sampling schedules are just data frames and can be merged easily
samples <- rbind(modern_samples, ancient_samples)

# run a simulation using the msprime back end from a compiled slendr model object
ts <- msprime(model, sequence_length = 1e5, recombination_rate = 0, samples = samples)

# simulated tree-sequence object can be saved to a file using ts_write()...
ts_file <- normalizePath(tempfile(fileext = ".trees"), winslash = "/", mustWork = FALSE)
ts_write(ts, ts_file)
# ... and, at a later point, loaded by ts_read()
ts <- ts_read(ts_file, model)

summary(ts)

Description

Generate the overlap of two slendr objects

Usage

overlap(x, y, name = NULL)

Arguments

Value

Object of the class slendr_region which encodes a standard spatial object of the class sf with several additional attributes (most importantly a corresponding slendr_map object, if applicable).

Description

Plots objects of the three slendr spatial classes (slendr_map, slendr_region, and slendr_pop).

Usage

plot_map(
  ...,
  time = NULL,
  gene_flow = FALSE,
  graticules = "original",
  intersect = TRUE,
  show_map = TRUE,
  title = NULL,
  interpolated_maps = NULL
)

Arguments

Value

A ggplot2 object with the visualized slendr map

Description

Plot demographic history encoded in a slendr model

Usage

plot_model(
  model,
  sizes = TRUE,
  proportions = FALSE,
  gene_flow = TRUE,
  log = FALSE,
  order = NULL,
  file = NULL,
  samples = NULL,
  ...
)

Arguments

Value

A ggplot2 object with the visualized slendr model

Examples

init_env()

# load an example model with an already simulated tree sequence
path <- system.file("extdata/models/introgression", package = "slendr")
model <- read_model(path)

plot_model(model, sizes = FALSE, log = TRUE)

Description

Defines the parameters of a population (non-spatial and spatial).

Usage

population(
  name,
  time,
  N,
  parent = NULL,
  map = FALSE,
  center = NULL,
  radius = NULL,
  polygon = NULL,
  remove = NULL,
  intersect = TRUE,
  competition = NA,
  mating = NA,
  dispersal = NA,
  dispersal_fun = NULL,
  aquatic = FALSE
)

Arguments

Details

There are four ways to specify a spatial boundary: i) circular range specified using a center coordinate and a radius, ii) polygon specified as a list of two-dimensional vector coordinates, iii) polygon as in ii), but defined (and named) using the region function, iv) with just a world map specified (circular or polygon range parameters set to the default NULL value), the population will be allowed to occupy the entire landscape.

Note that because slendr models have to accomodate both SLiM and msprime back ends, population sizes and split times are rounded to the nearest integer value.

Value

Object of the class slendr_pop, which contains population parameters such as name, time of appearance in the simulation, parent population (if any), and its spatial parameters such as map and spatial boundary.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

All spatial objects in the slendr package are internally represented as Simple Features (sf) objects. This fact is hidden in most circumstances this, as the goal of the slendr package is to provide functionality at a much higher level (population boundaries, geographic regions, instead of individual polygons and other "low-level" geometric objects), without the users having to worry about low-level details involved in handling spatial geometries. However, the full sf object representation can be always printed by calling x[].

Usage

## S3 method for class 'slendr_pop'
print(x, ...)

## S3 method for class 'slendr_region'
print(x, ...)

## S3 method for class 'slendr_map'
print(x, ...)

## S3 method for class 'slendr_model'
print(x, ...)

Arguments

Value

No return value, used only for printing

Description

Print tskit's summary table of the Python tree-sequence object

Usage

## S3 method for class 'slendr_ts'
print(x, ...)

Arguments

Value

No return value, simply prints the tskit summary table to the terminal

Description

Reads all configuration tables and other model data from a location where it was previously compiled to by the compile function.

Usage

read_model(path)

Arguments

Value

Compiled slendr_model model object which encapsulates all information about the specified model (which populations are involved, when and how much gene flow should occur, what is the spatial resolution of a map, and what spatial dispersal and mating parameters should be used in a SLiM simulation, if applicable)

Examples

init_env()

# load an example model with an already simulated tree sequence
path <- system.file("extdata/models/introgression", package = "slendr")
model <- read_model(path)

plot_model(model, sizes = FALSE, log = TRUE)

Description

Creates a geographic region (a polygon) on a given map and gives it a name. This can be used to define objects which can be reused in multiple places in a slendr script (such as region arguments of population) without having to repeatedly define polygon coordinates.

Usage

region(name = NULL, map = NULL, center = NULL, radius = NULL, polygon = NULL)

Arguments

Value

Object of the class slendr_region which encodes a standard spatial object of the class sf with several additional attributes (most importantly a corresponding slendr_map object, if applicable).

Examples

# create a blank abstract world 1000x1000 distance units in size
blank_map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# it is possible to construct custom landscapes (islands, corridors, etc.)
island1 <- region("island1", polygon = list(c(10, 30), c(50, 30), c(40, 50), c(0, 40)))
island2 <- region("island2", polygon = list(c(60, 60), c(80, 40), c(100, 60), c(80, 80)))
island3 <- region("island3", center = c(20, 80), radius = 10)
archipelago <- island1 %>% join(island2) %>% join(island3)

custom_map <- world(xrange = c(1, 100), c(1, 100), landscape = archipelago)

# real Earth landscapes can be defined using freely-available Natural Earth
# project data and with the possibility to specify an appropriate Coordinate
# Reference System, such as this example of a map of Europe

real_map <- world(xrange = c(-15, 40), yrange = c(30, 60), crs = "EPSG:3035")

Description

Converts between coordinates on a compiled raster map (i.e. pixel units) and different Geographic Coordinate Systems (CRS).

Usage

reproject(
  from,
  to,
  x = NULL,
  y = NULL,
  coords = NULL,
  model = NULL,
  add = FALSE,
  input_prefix = "",
  output_prefix = "new"
)

Arguments

Value

Data.frame with converted two-dimensional coordinates given as input

Examples

lon_lat_df <- data.frame(x = c(30, 0, 15), y = c(60, 40, 10))

reproject(
  from = "epsg:4326",
  to = "epsg:3035",
  coords = lon_lat_df,
  add = TRUE # add converted [lon,lat] coordinates as a new column
)

Description

Resizes the population starting from the current value of N individuals to the specified value

Usage

resize(pop, N, how, time, end = NULL)

Arguments

Details

In the case of exponential size change, if the final N is larger than the current size, the population will be exponentially growing over the specified time period until it reaches N individuals. If N is smaller, the population will shrink exponentially.

Note that because slendr models have to accomodate both SLiM and msprime back ends, population sizes and split times are rounded to the nearest integer value.

Value

Object of the class slendr_pop, which contains population parameters such as name, time of appearance in the simulation, parent population (if any), and its spatial parameters such as map and spatial boundary.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

Schedule sampling events at specified times and, optionally, a given set of locations on a landscape

Usage

schedule_sampling(model, times, ..., locations = NULL, strict = FALSE)

Arguments

Details

If both times and locations are given, the the sampling will be scheduled on each specified location in each given time-point. Note that for the time-being, in the interest of simplicity, no sanity checks are performed on the locations given except the restriction that the sampling points must fall within the bounding box around the simulated world map. Other than that, slendr will simply instruct its SLiM backend script to sample individuals as close to the sampling points given as possible, regardless of whether those points lie within a population spatial boundary at that particular moment of time.

Value

Data frame with three columns: time of sampling, population to sample from, how many individuals to sample

Examples

init_env()

# load an example model with an already simulated tree sequence
path <- system.file("extdata/models/introgression", package = "slendr")
model <- read_model(path)

# afr and eur objects would normally be created before slendr model compilation,
# but here we take them out of the model object already compiled for this
# example (in a standard slendr simulation pipeline, this wouldn't be necessary)
afr <- model$populations[["AFR"]]
eur <- model$populations[["EUR"]]

# schedule the recording of 10 African and 100 European individuals from a
# given model at 20 ky, 10 ky, 5ky ago and at present-day (time 0)
schedule <- schedule_sampling(
  model, times = c(20000, 10000, 5000, 0),
  list(afr, 10), list(eur, 100)
)

# the result of `schedule_sampling` is a simple data frame (note that the locations
# of sampling locations have `NA` values because the model is non-spatial)
schedule

Description

Changes either the competition interactive distance, mating choice distance, or the dispersal of offspring from its parent

Usage

set_dispersal(
  pop,
  time,
  competition = NA,
  mating = NA,
  dispersal = NA,
  dispersal_fun = NULL
)

Arguments

Value

Object of the class slendr_pop, which contains population parameters such as name, time of appearance in the simulation, parent population (if any), and its spatial parameters such as map and spatial boundary.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

This function allows a more manual control of spatial map changes in addition to the expand and move functions

Usage

set_range(
  pop,
  time,
  center = NULL,
  radius = NULL,
  polygon = NULL,
  lock = FALSE
)

Arguments

Value

Object of the class slendr_pop, which contains population parameters such as name, time of appearance in the simulation, parent population (if any), and its spatial parameters such as map and spatial boundary.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

This function will automatically download a Python miniconda distribution dedicated to an R-Python interface. It will also create a slendr-specific Python environment with all the required Python dependencies.

Usage

setup_env(quiet = FALSE, agree = FALSE, pip = FALSE)

Arguments

Value

No return value, called for side effects

Description

Shrinks the spatial population range by a specified distance in a given time-window

Usage

shrink_range(
  pop,
  by,
  end,
  start,
  overlap = 0.8,
  snapshots = NULL,
  lock = FALSE,
  verbose = TRUE
)

Arguments

Details

Note that because slendr models have to accomodate both SLiM and msprime back ends, population sizes and split times are rounded to the nearest integer value.

Value

Object of the class slendr_pop, which contains population parameters such as name, time of appearance in the simulation, parent population (if any), and its spatial parameters such as map and spatial boundary.

Examples

# spatial definitions -----------------------------------------------------

# create a blank abstract world 1000x1000 distance units in size
map <- world(xrange = c(0, 1000), yrange = c(0, 1000), landscape = "blank")

# create a circular population with the center of a population boundary at
# [200, 800] and a radius of 100 distance units, 1000 individuals at time 1
# occupying a map just specified
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100)

# printing a population object to a console shows a brief summary
pop1

# create another population occupying a polygon range, splitting from pop1
# at a given time point (note that specifying a map is not necessary because
# it is "inherited" from the parent)
pop2 <- population("pop2", N = 100, time = 50, parent = pop1,
                        polygon = list(c(100, 100), c(320, 30), c(500, 200),
                                  c(500, 400), c(300, 450), c(100, 400)))

pop3 <- population("pop3", N = 200, time = 80, parent = pop2,
                   center = c(800, 800), radius = 200)

# move "pop1" to another location along a specified trajectory and saved the
# resulting object to the same variable (the number of intermediate spatial
# snapshots can be also determined automatically by leaving out the
# `snapshots = ` argument)
pop1_moved <- move(pop1, start = 100, end = 200, snapshots = 6,
                   trajectory = list(c(600, 820), c(800, 400), c(800, 150)))
pop1_moved

# many slendr functions are pipe-friendly, making it possible to construct
# pipelines which construct entire history of a population
pop1 <- population("pop1", N = 1000, time = 1,
                   map = map, center = c(200, 800), radius = 100) %>%
  move(start = 100, end = 200, snapshots = 6,
       trajectory = list(c(400, 800), c(600, 700), c(800, 400), c(800, 150))) %>%
  set_range(time = 300, polygon = list(
    c(400, 0), c(1000, 0), c(1000, 600), c(900, 400), c(800, 250),
    c(600, 100), c(500, 50))
  )

# population ranges can expand by a given distance in all directions
pop2 <- expand_range(pop2, by = 200, start = 50, end = 150, snapshots = 3)

# we can check the positions of all populations interactively by plotting their
# ranges together on a single map
plot_map(pop1, pop2, pop3)

# gene flow events --------------------------------------------------------

# individual gene flow events can be saved to a list
gf <- list(
  gene_flow(from = pop1, to = pop3, start = 150, end = 200, rate = 0.15),
  gene_flow(from = pop1, to = pop2, start = 300, end = 330, rate = 0.25)
)

# compilation -------------------------------------------------------------

# compile model components in a serialized form to dist, returning a single
# slendr model object (in practice, the resolution should be smaller)
model <- compile_model(
  populations = list(pop1, pop2, pop3), generation_time = 1,
  resolution = 100, simulation_length = 500,
  competition = 5, mating = 5, dispersal = 1
)

Description

This function will execute a SLiM script generated by the compile function during the compilation of a slendr demographic model.

Usage

slim(
  model,
  sequence_length,
  recombination_rate,
  samples = NULL,
  ts = TRUE,
  path = NULL,
  random_seed = NULL,
  method = c("batch", "gui"),
  verbose = FALSE,
  run = TRUE,
  slim_path = NULL,
  burnin = 0,
  max_attempts = 1,
  spatial = !is.null(model$world),
  coalescent_only = TRUE,
  locations = NULL
)

Arguments

Details

The arguments sequence_length and recombination_rate can be omitted for slendr models utilizing customized initialization of genomic architecture. In such cases, users may either provide hard-coded values directly through SLiM's initializeGenomicElement() and initializeRecombinationRate() functions or utilize slendr's templating functionality provided by its substitute() function. When ts = TRUE, the returning value of this function depends on whether or not the path argument was set. If the user did provide the path where output files should be saved, the path is returned (invisibly). This is mostly intended to support simulations of customized user models. If path is not set by the user, it is assumed that a tree-sequence object is desired as a sole return value of the function (when ts = TRUE) and so it is automatically loaded when simulation finishes, or (when ts = FALSE) that only customized files are to be produced by the simulation, in which the user will be loading such files by themselves (and only the path is needed).

Value

A tree-sequence object loaded via Python-R reticulate interface function ts_read (internally represented by the Python object tskit.trees.TreeSequence). If the path argument was set, specifying the directory where results should be saved, the function will return this path as a single-element character vector.

Examples

init_env()

# load an example model
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# afr and eur objects would normally be created before slendr model compilation,
# but here we take them out of the model object already compiled for this
# example (in a standard slendr simulation pipeline, this wouldn't be necessary)
afr <- model$populations[["AFR"]]
eur <- model$populations[["EUR"]]
chimp <- model$populations[["CH"]]

# schedule the sampling of a couple of ancient and present-day individuals
# given model at 20 ky, 10 ky, 5ky ago and at present-day (time 0)
modern_samples <- schedule_sampling(model, times = 0, list(afr, 5), list(eur, 5), list(chimp, 1))
ancient_samples <- schedule_sampling(model, times = c(30000, 20000, 10000), list(eur, 1))

# sampling schedules are just data frames and can be merged easily
samples <- rbind(modern_samples, ancient_samples)

# run a simulation using the SLiM back end from a compiled slendr model object and return
# a tree-sequence object as a result
ts <- slim(model, sequence_length = 1e5, recombination_rate = 0, samples = samples)

# simulated tree-sequence object can be saved to a file using ts_write()...
ts_file <- normalizePath(tempfile(fileext = ".trees"), winslash = "/", mustWork = FALSE)
ts_write(ts, ts_file)
# ... and, at a later point, loaded by ts_read()
ts <- ts_read(ts_file, model)

ts

Description

Substitute values of templated {{parameters}} in a given SLiM extension template

Usage

substitute_values(template, ...)

Arguments

Details

If a file or a multi-line string given as template contains parameters specified as {{param}} where "param" can be arbitrary variable name, this function substitutes each templated {{parameter}} for a given values. Such modified template is then used to extend a built-in slendr SLiM script, allowing for a customization of its default behavior (most commonly replacing its assumption of neutrality for non-neutral scenarios, such as simulations of natural selection).

Value

Path to a file with a saved extension script containing all substituted values

Description

Generate the difference between two slendr objects

Usage

subtract(x, y, name = NULL)

Arguments

Value

Object of the class slendr_region which encodes a standard spatial object of the class sf with several additional attributes (most importantly a corresponding slendr_map object, if applicable).

Description

Summarise the contents of a ts_nodes result

Usage

## S3 method for class 'slendr_nodes'
summary(object, ...)

Arguments

Value

Used for its output to the terminal

Description

This function computes the AFS with respect to the given set of individuals or nodes.

Usage

ts_afs(
  ts,
  sample_sets = NULL,
  mode = c("site", "branch", "node"),
  windows = NULL,
  span_normalise = FALSE,
  polarised = TRUE
)

Arguments

Details

For more information on the format of the result and dimensions, in particular the interpretation of the first and the last element of the AFS (when complete = TRUE), please see the tskit manual at https://tskit.dev/tskit/docs/stable/python-api.html and the example section dedicated to AFS at https://tskit.dev/tutorials/analysing_tree_sequences.html#allele-frequency-spectra.

Value

Allele frequency spectrum values for the given sample set. Note that the contents of the first and last elements of the AFS might surprise you. Read the links in the description for more detail on how tskit handles things.

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model) %>% ts_mutate(mutation_rate = 1e-8, random_seed = 42)

samples <- ts_samples(ts) %>% .[.$pop %in% c("AFR", "EUR"), ]

# compute AFS for the given set of individuals
ts_afs(ts, sample_sets = list(samples$name))

Description

Extract (spatio-)temporal ancestral history for given nodes/individuals

Usage

ts_ancestors(ts, x, verbose = FALSE, complete = TRUE)

Arguments

Value

A table of ancestral nodes of a given tree-sequence node all the way up to the root of the tree sequence

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model)

# find the complete ancestry information for a given individual
ts_ancestors(ts, "EUR_1", verbose = TRUE)

Description

Check that all trees in the tree sequence are fully coalesced

Usage

ts_coalesced(ts, return_failed = FALSE)

Arguments

Value

TRUE or FALSE value if return_failed = FALSE, otherwise a vector of (tskit Python 0-based) indices of trees which failed the coalescence test

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model)

ts_coalesced(ts) # is the tree sequence fully coalesced? (TRUE or FALSE)

# returns a vector of tree sequence segments which are not coalesced
not_coalesced <- ts_coalesced(ts, return_failed = TRUE)

Description

Extract all descendants of a given tree-sequence node

Usage

ts_descendants(ts, x, verbose = FALSE, complete = TRUE)

Arguments

Value

A table of descendant nodes of a given tree-sequence node all the way down to the leaves of the tree sequence

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model)

# find the complete descendancy information for a given individual
ts_descendants(ts, x = 62, verbose = TRUE)

Description

Calculate pairwise divergence between sets of individuals

Usage

ts_divergence(
  ts,
  sample_sets,
  mode = c("site", "branch", "node"),
  windows = NULL,
  span_normalise = TRUE
)

Arguments

Value

For each pairwise calculation, either a single divergence value or a vector of divergence values (one for each window)

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model) %>% ts_mutate(mutation_rate = 1e-8, random_seed = 42)

# collect sampled individuals from all populations in a list
sample_sets <- ts_samples(ts) %>%
  split(., .$pop) %>%
  lapply(function(pop) pop$name)

# compute the divergence between individuals from each sample set (list of
# individual names generated in the previous step)
ts_divergence(ts, sample_sets) %>% .[order(.$divergence), ]

Description

Calculate diversity in given sets of individuals

Usage

ts_diversity(
  ts,
  sample_sets,
  mode = c("site", "branch", "node"),
  windows = NULL,
  span_normalise = TRUE
)

Arguments

Value

For each set of individuals either a single diversity value or a vector of diversity values (one for each window)

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model) %>% ts_mutate(mutation_rate = 1e-8, random_seed = 42)

# collect sampled individuals from all populations in a list
sample_sets <- ts_samples(ts) %>%
  split(., .$pop) %>%
  lapply(function(pop) pop$name)

# compute diversity in each population based on sample sets extracted
# in the previous step
ts_diversity(ts, sample_sets) %>% .[order(.$diversity), ]

Description

This function first obtains an SVG representation of the tree by calling the draw_svg method of tskit and renders it as a bitmap image in R. All of the many optional keyword arguments of the draw_svg method can be provided and will be automatically passed to the method behind the scenes.

Usage

ts_draw(
  x,
  width = 1000,
  height = 1000,
  labels = FALSE,
  sampled_only = TRUE,
  title = NULL,
  ...
)

Arguments

Value

No return value, called for side effects

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model)

# extract the first tree in the tree sequence and draw it
tree <- ts_tree(ts, i = 1)

# ts_draw accepts various optional arguments of tskit.Tree.draw_svg
ts_draw(tree, time_scale = "rank")

Description

Extract spatio-temporal edge annotation table from a given tree or tree sequence

Usage

ts_edges(x)

Arguments

Value

Data frame of the sf type containing the times of nodes and start-end coordinates of edges across space

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model)

# extract an annotated table with (spatio-)temporal edge information
ts_edges(ts)

Description

EIGENSTRAT data produced by this function can be used by the admixr R package (https://bodkan.net/admixr/).

Usage

ts_eigenstrat(ts, prefix, chrom = "chr1", outgroup = NULL)

Arguments

Details

In case an outgroup was not formally specified in a slendr model which generated the tree sequence data, it is possible to artificially create an outgroup sample with the name specified by the outgroup argument, which will carry all ancestral alleles (i.e. value "2" in a geno file for each position in a snp file).

Value

Object of the class EIGENSTRAT created by the admixr package

Description

These functions present an R interface to the corresponding f-statistics methods in tskit.

Usage

ts_f2(
  ts,
  A,
  B,
  mode = c("site", "branch", "node"),
  span_normalise = TRUE,
  windows = NULL
)

ts_f3(
  ts,
  A,
  B,
  C,
  mode = c("site", "branch", "node"),
  span_normalise = TRUE,
  windows = NULL
)

ts_f4(
  ts,
  W,
  X,
  Y,
  Z,
  mode = c("site", "branch", "node"),
  span_normalise = TRUE,
  windows = NULL
)

ts_f4ratio(
  ts,
  X,
  A,
  B,
  C,
  O,
  mode = c("site", "branch"),
  span_normalise = TRUE
)

Arguments

Details

Note that the order of populations f3 statistic implemented in tskit (https://tskit.dev/tskit/docs/stable/python-api.html#tskit.TreeSequence.f3) is different from what you might expect from ADMIXTOOLS, as defined in Patterson 2012 (see doi:10.1534/genetics.112.145037 under heading "The three-population test and introduction of f-statistics", as well as ADMIXTOOLS documentation at https://github.com/DReichLab/AdmixTools/blob/master/README.3PopTest#L5). Specifically, the widely used notation introduced by Patterson assumes the population triplet as f3(C; A, B), with C being the "focal" sample (i.e., either the outgroup or a sample tested for admixture). In contrast, tskit implements f3(A; B, C), with the "focal sample" being A.

Although this is likely to confuse many ADMIXTOOLS users, slendr does not have much choice in this, because its ts_*() functions are designed to be broadly compatible with raw tskit methods.

Value

Data frame with statistics calculated for the given sets of individuals

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk and add mutations to it
ts <- ts_read(slendr_ts, model) %>% ts_mutate(mutation_rate = 1e-8, random_seed = 42)

# calculate f2 for two individuals in a previously loaded tree sequence
ts_f2(ts, A = "AFR_1", B = "EUR_1")

# calculate f2 for two sets of individuals
ts_f2(ts, A = c("AFR_1", "AFR_2"), B = c("EUR_1", "EUR_3"))

# calculate f3 for two individuals in a previously loaded tree sequence
ts_f3(ts, A = "EUR_1", B = "AFR_1", C = "NEA_1")

# calculate f3 for two sets of individuals
ts_f3(ts, A = c("AFR_1", "AFR_2", "EUR_1", "EUR_2"),
          B = c("NEA_1", "NEA_2"),
          C = "CH_1")

# calculate f4 for single individuals
ts_f4(ts, W = "EUR_1", X = "AFR_1", Y = "NEA_1", Z = "CH_1")

# calculate f4 for sets of individuals
ts_f4(ts, W = c("EUR_1", "EUR_2"),
          X = c("AFR_1", "AFR_2"),
          Y = "NEA_1",
          Z = "CH_1")

# calculate f4-ratio for a given set of target individuals X
ts_f4ratio(ts, X = c("EUR_1", "EUR_2", "EUR_4", "EUR_5"),
               A = "NEA_1", B = "NEA_2", C = "AFR_1", O = "CH_1")

Description

For a discussion on the difference between "site", "branch", and "node" options of the mode argument, please see the tskit documentation at https://tskit.dev/tskit/docs/stable/stats.html#sec-stats-mode.

Usage

ts_fst(
  ts,
  sample_sets,
  mode = c("site", "branch", "node"),
  windows = NULL,
  span_normalise = TRUE
)

Arguments

Value

For each pairwise calculation, either a single Fst value or a vector of Fst values (one for each window)

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model) %>% ts_mutate(mutation_rate = 1e-8, random_seed = 42)

# compute F_st between two sets of individuals in a given tree sequence ts
ts_fst(ts, sample_sets = list(afr = c("AFR_1", "AFR_2", "AFR_3"),
                              eur = c("EUR_1", "EUR_2")))

Description

Extract genotype table from the tree sequence

Usage

ts_genotypes(ts)

Arguments

Value

Data frame object of the class tibble containing genotypes of simulated individuals in columns

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk, recapitate it, simplify it, and mutate it
ts <- ts_read(slendr_ts, model) %>%
  ts_recapitate(Ne = 10000, recombination_rate = 1e-8) %>%
  ts_simplify() %>%
  ts_mutate(mutation_rate = 1e-8)

# extract the genotype matrix (this could take  a long time consume lots
# of memory!)
gts <- ts_genotypes(ts)

Description

This function iterates over a tree sequence and returns IBD tracts between pairs of individuals or nodes

Usage

ts_ibd(
  ts,
  coordinates = FALSE,
  within = NULL,
  between = NULL,
  squash = FALSE,
  minimum_length = NULL,
  maximum_time = NULL,
  sf = TRUE
)

Arguments

Details

This function is considered experimental. For full control over IBD segment detection in tree-sequence data, users can (and perhaps, for the time being, should) rely on the tskit method ibd_segments (see https://tskit.dev/tskit/docs/stable/python-api.html#tskit.TreeSequence.ibd_segments).

Iternally, this function leverages the tskit TreeSequence method ibd_segments. However, note that the ts_ibd function always returns a data frame of IBD tracts, it does not provide an option to iterate over individual IBD segments as shown in the official tskit documentation at https://tskit.dev/tskit/docs/stable/ibd.html. In general, R handles heavy iteration poorly, and this function does not attempt to serve as a full wrapper to ibd_segments.

Unfortunately, the distinction between "squashed IBD" (what many would consider to be the expected definition of IBD) and tskit’s IBD which is defined via distinct genealogical paths (see https://github.com/tskit-dev/tskit/issues/2459 for a discussion of the topic), makes the meaning of the filtering parameter of the ibd_segments() method of tskit minimum_length somewhat unintuitive. As of this moment, this function argument filters on IBD segments on the tskit level, not the level of the squashed IBD segments!

Value

A data frame with IBD results (either coordinates of each IBD segment shared by a pair of nodes, or summary statistics about the total IBD sharing for that pair)

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model)

# find IBD segments between specified Neanderthals and Europeans
ts_ibd(
  ts,
  coordinates = TRUE,
  between = list(c("NEA_1", "NEA_2"), c("EUR_1", "EUR_2")),
  minimum_length = 40000
)

Description

Deprecated function. Please use ts_read instead.

Usage

ts_load(file, model = NULL)

Arguments

Description

Extract list with tree sequence metadata saved by SLiM

Usage

ts_metadata(ts)

Arguments

Value

List of metadata fields extracted from the tree-sequence object

Examples

init_env()

# load an example model with an already simulated tree sequence
slendr_ts <- system.file("extdata/models/introgression_slim.trees", package = "slendr")
model <- read_model(path = system.file("extdata/models/introgression", package = "slendr"))

# load the tree-sequence object from disk
ts <- ts_read(slendr_ts, model)

# extract the list of metadata information from the tree sequence
ts_metadata(ts)

Description

Add mutations to the given tree sequence

Usage

ts_mutate(
  ts,
  mutation_rate,
  random_seed = NULL,
  keep_existing = TRUE,
  mutation_model = NULL
)

Arguments

Value

Tree-sequence object of the class slendr_ts, which serves as an interface point for the Python module tskit using slendr functions with the ts_ prefix.

See Also

ts_nodes for extracting useful information about individuals, nodes, coalescent times and geospatial locations of nodes on a map