In thibautjombart/quicksim: Quick and simple outbreak simulator

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  fig.path = "figs/",
  echo = TRUE
)

Installing the package

To install the development version of the package, type:

devtools::install_github("thibautjombart/quicksim")

Note that this requires the package devtools installed.

What does it do?

As it currently stands, the package does not implement the process generating the transmission tree, but it does provide functions to create cases given known transmission events.

The main features of the package include:

new_case: a function creating a new case, either from a prior infector, or imported
new_location: a function creating a new case location, either from a prior infector, or imported; used in new_case
new_dna: a function creating the DNA sequence for a new case, either from a prior infector DNA sequence, or imported; used in new_case
dist_dna: a function to compute pairwise genetic distances from 'case' data

Demo

A basic example

We start from the following transmission tree:

tree <- matrix(c(1, 2, 2, 3, 3, 4, 2, 5, 5, 6), ncol = 2, byrow = TRUE)
colnames(tree) <- c("from", "to")
tree

## dates of infection for cases 1-6
infection_dates <- c(0, 2, 5, 6, 7, 10)

library(igraph)
plot(graph.data.frame(tree), vertex.color = "#e0ccff")

We load the package, examine the default settings, make custom config, and generate an index case:

library(quicksim)

## see default config
new_config()

## use custom config
conf <- new_config(genome_length = 1e4,
                   mutation_rate = 1e-4,
                   separation_lineages = 10)

## create index case
index <- new_case(config = conf)
index

We create new cases according to the transmission tree:

set.seed(1)

cases <- vector(mode = "list", length = 6)
cases[[1]] <- index

for (i in seq_len(nrow(tree))) {
  infector <- tree[i, 1]
  infectee <- tree[i, 2]
  cases[[infectee]] <- new_case(cases[[infector]],
                                date = infection_dates[infectee],
                                config = conf)
}

cases

Some further analyses

library(ape)

dna <- lapply(cases, function(e) e$dna)
locations <- matrix(unlist(lapply(cases, function(e) e$location)),
                    ncol = 2, byrow = TRUE)

## genetic distances
D_dna <- dist_dna(dna)
D_dna

plot(nj(D_dna), "unrooted", main = "Unrooted NJ tree, Hamming distances",
     show.tip.label = FALSE)

tiplabels(frame = "circ", col = "#f9ecf2", bg = "#602040", cex = 1.5)

## spatial distances
D_geo <- dist(locations)
D_geo
plot(locations, pch = paste(1:6), xlab = "x", ylab = "y", cex = 2)

Using custom spatial kernel

Customised spatial kernels can be specified as part of the config through the argument spatial_kernel. The behaviour of new_location is as follows:

by default, a Normal kernel with standard deviation 1 is used
if provided, sd_spatial is used for the standard deviation of the Normal kernel
if provided, the alternative spatial kernel spatial_kernel is used instead, in which case sd_spatial is ignored.

If provided, the argument spatial_kernel should be a function of x, a vector of spatial coordinates, and return a similarly sized vector of new coordinates.

We can repeat the previous example with a new spatial kernel:

## new custom config
custom_kernel <- function(x) {
  out <- x + runif(length(x), min = 1, max = 3)
  out
}

new_conf <- new_config(genome_length = 1e4,
                       mutation_rate = 1e-4,
                       separation_lineages = 10,
                       spatial_kernel = custom_kernel)

set.seed(1)

for (i in seq_len(nrow(tree))) {
  infector <- tree[i, 1]
  infectee <- tree[i, 2]
  cases[[infectee]] <- new_case(cases[[infector]],
                                date = infection_dates[infectee],
                                config = new_conf)
}


## look at the new spatial distribution
locations <- matrix(unlist(lapply(cases, function(e) e$location)),
                    ncol = 2, byrow = TRUE)

plot(locations, pch = paste(1:6), xlab = "x", ylab = "y", cex = 2)

Resources

Vignettes

Vignettes are not currently available for this package.

Websites

The following websites are available:

The quicksim project on github, useful for developers, contributors, and users wanting to post issues, bug reports and feature requests:
http://github.com/thibautjombart/quicksim

Getting help online

Bug reports and feature requests should be posted on github using the issue system. All other questions should be posted on the RECON forum:
http://www.repidemicsconsortium.org/forum/