README.md

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:

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
#>      from to
#> [1,]    1  2
#> [2,]    2  3
#> [3,]    3  4
#> [4,]    2  5
#> [5,]    5  6

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

library(igraph)
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union
plot(graph.data.frame(tree), vertex.color = "#e0ccff")

plot of chunk tree

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

library(quicksim)

## see default config
new_config()
#> $x_min
#> [1] 0
#> 
#> $x_max
#> [1] 100
#> 
#> $y_min
#> [1] 0
#> 
#> $y_max
#> [1] 100
#> 
#> $spatial_kernel
#> function (x) 
#> stats::rnorm(length(x), mean = x, sd = sd)
#> <environment: 0x31b8a90>
#> 
#> $sd_spatial
#> [1] 1
#> 
#> $genome_length
#> [1] 30000
#> 
#> $mutation_rate
#> [1] 1e-05
#> 
#> $separation_lineages
#> [1] 365

## 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
#> <case object>
#> 
#> /// date of infection ($date): 
#> [1] 0
#> 
#> /// place of infection ($location): 
#> [1] 96.45478 32.31864
#> 
#> /// DNA sequence mutations ($dna): 
#> [1] 9703 9129 9428 6714 6923 1048 6664 5071 1029

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
#> [[1]]
#> <case object>
#> 
#> /// date of infection ($date): 
#> [1] 0
#> 
#> /// place of infection ($location): 
#> [1] 96.45478 32.31864
#> 
#> /// DNA sequence mutations ($dna): 
#> [1] 9703 9129 9428 6714 6923 1048 6664 5071 1029
#> 
#> [[2]]
#> <case object>
#> 
#> /// date of infection ($date): 
#> [1] 2
#> 
#> /// place of infection ($location): 
#> [1] 96.63842 31.48301
#> 
#> /// DNA sequence mutations ($dna): 
#>  [1] 9703 9129 9428 6714 6923 1048 6664 5071 1029 3722
#> 
#> [[3]]
#> <case object>
#> 
#> /// date of infection ($date): 
#> [1] 5
#> 
#> /// place of infection ($location): 
#> [1] 96.34370 31.47725
#> 
#> /// DNA sequence mutations ($dna): 
#>  [1] 9703 9129 9428 6714 6923 1048 6664 5071 1029 3722 6608 6292  619 2061
#> [15] 1766 6871
#> 
#> [[4]]
#> <case object>
#> 
#> /// date of infection ($date): 
#> [1] 6
#> 
#> /// place of infection ($location): 
#> [1] 96.73355 30.85601
#> 
#> /// DNA sequence mutations ($dna): 
#>  [1] 9703 9129 9428 6714 6923 1048 6664 5071 1029 3722 6608 6292  619 2061
#> [15] 1766 6871 3801 7775 9347 2122
#> 
#> [[5]]
#> <case object>
#> 
#> /// date of infection ($date): 
#> [1] 7
#> 
#> /// place of infection ($location): 
#> [1] 97.76335 31.43808
#> 
#> /// DNA sequence mutations ($dna): 
#>  [1] 9703 9129 9428 6714 6923 1048 6664 5071 1029 3722 3824
#> 
#> [[6]]
#> <case object>
#> 
#> /// date of infection ($date): 
#> [1] 10
#> 
#> /// place of infection ($location): 
#> [1] 98.58458 32.03198
#> 
#> /// DNA sequence mutations ($dna): 
#>  [1] 9703 9129 9428 6714 6923 1048 6664 5071 1029 3722 3824 1863 8274 6685

Some further analyses

library(ape)
#> 
#> Attaching package: 'ape'
#> The following objects are masked from 'package:igraph':
#> 
#>     edges, mst, ring

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
#>    1  2  3  4  5
#> 2  1            
#> 3  7  6         
#> 4 11 10  4      
#> 5  2  1  7 11   
#> 6  5  4 10 14  3

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)

plot of chunk analyses


## spatial distances
D_geo <- dist(locations)
D_geo
#>           1         2         3         4         5
#> 2 0.8555700                                        
#> 3 0.8486961 0.2947769                              
#> 4 1.4889646 0.6341822 0.7334287                    
#> 5 1.5772623 1.1258280 1.4201915 1.1829265          
#> 6 2.1490005 2.0220963 2.3085148 2.1929952 1.0134708
plot(locations, pch = paste(1:6), xlab = "x", ylab = "y", cex = 2)

plot of chunk analyses

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:

  1. by default, a Normal kernel with standard deviation 1 is used
  2. if provided, sd_spatial is used for the standard deviation of the Normal kernel
  3. 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)

plot of chunk custom_kernel

Resources

Vignettes

Vignettes are not currently available for this package.

Websites

The following websites are available:

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/

A quick overview

The following worked example provides a brief overview of the package's functionalities.

Contributors (by alphabetic order):

See details of contributions on: https://github.com/thibautjombart/quicksim/graphs/contributors

Contributions are welcome via pull requests.

Please note that this project is released with a Contributor Code of Conduct. By participating in this project you agree to abide by its terms.

Maintainer: Thibaut Jombart (thibautjombart@gmail.com)



thibautjombart/quicksim documentation built on May 5, 2019, 2:42 a.m.