Analyzing non-slendr tree sequences

Non-spatial SLiM tree sequences

Consider the following SLiM script, which creates a couple of populations (with different N_e) splitting from an ancestral population p1 (lets save it to /tmp/nonspatial.slim):

After we run this script in SLiM, we can use slendr to load the output tree sequence (saved to /tmp/nonspatial-slim.trees), simplify it, and overlay mutations on it using the standard functionality originally developed for slendr tree sequences. Note that this is the same command we would use for loading slendr tree sequences, except we direct the ts_read() function straight to the tree-sequence output file rather than using the ts_read(<model-object>) format used when working with standard slendr simulations. This way, slendr

ts <- ts_read(nonspatial_trees_file) %>%
  ts_simplify() %>%
  ts_mutate(mutation_rate = 1e-7, random_seed = SEED)

We can extract information about individual’s names, nodes, population assignments, etc. just as with any slendr tree sequence with the function ts_nodes(). As with standard slendr models, this function loads the “raw” node and individual tree-sequences tables, performs a couple of join operations, and presents the whole thing in a nice unified form for interactive data analysis (which can also include spatial information—see below):

data <- ts_nodes(ts) %>% dplyr::filter(sampled)
data

Moving on to tskit statistics, we can use the data table above to extract a list of nodes belonging to each population (this is what various tskit tree-sequence statistics operate on, and slendr follows that design). Here we are computing the nucleotide diversity in each of the four populations using the ts_diversity() function, by first creating a list of lists with node IDs (i.e. chromosomes) of individuals in each population:

sample_sets <- split(data$node_id, data$pop)

# compute nucleotide diversity in each population
# (any other ts_*() tskit R interface function should work)
ts_diversity(ts, sample_sets)

Just as with slendr tree sequences (as demonstrated in our paper) we can get a individual trees too, extracted in the in the phylogenetic format provided by the ape R package. Here we first simplify the tree sequence even further to just 10 nodes to make things manageable:

samples <- sample(data$node_id, 10)
ts_small <- ts_simplify(ts, simplify_to = samples)

# extract the 42nd tree in the genealogy to an R 'phylo' format
tree <- ts_phylo(ts_small, 42 - 1)
tree

Once we have that R tree object, we can use packages like ggtree to visualize the tree (any other phylogenetic package would work too). Note that because nodes of ‘ape phylo’ trees must conform to a strict format (they must be labelled 1...N), we will extract the information about the node IDs in the tskit tree-sequence data to be able to plot them in the tree.

library(ggtree)

labels <- ts_nodes(tree) %>% select(node = phylo_id, tskit_id = node_id)

ggtree(tree, branch.length = "none") %<+% labels +
  geom_label(aes(label = tskit_id))

library(ape)
plot(tree, show.tip.label = FALSE)
nodelabels()
tiplabels()

msprime (non-slendr) tree sequences

The same as above applies also to msprime tree sequences (which is really not that surprising, given that it’s all tskit under the hood).

We can start with Python:

import msprime
ts = msprime.sim_ancestry(100)
ts.dump(<trees file>)

And then we can proceed with loading the msprime tree sequence into R and analyze it with the slendr functionality:

ts <- ts_read(msprime_trees_file)

ts_nodes(ts)

Spatial SLiM (non-slendr) tree sequences

Furthermore, the generalized interface also supports slendr’s spatial tree-sequence features, with all the bells and whistles.

For instance, lets take the following spatial SLiM script (modified from the SLiM manual) and execute it with SLiM the usual way:

We can then load and simplify the output tree sequence in just as we did above in this vignette (or anywhere in the slendr documentation):

ts <- ts_read(spatial_trees_file) %>% ts_simplify()

Finally, we can access the spatio-temporal data embedded in the output tree sequence in the standard slendr way (note the spatial sf column location with the POINT data type):

data <- ts_nodes(ts)
data

Because we get the tree sequence converted to the spatial sf data format, we can use standard geospatial packages to use any spatial data analysis methods that those packages provide.

To briefly demonstrate what this means, we can trivially plot the location of each recorded node:

ggplot() + geom_sf(data = data, aes(color = time), alpha = 0.5)

We can also collect spatio-temporal ancestry information of a particular node (i.e. the times and locations of all of its ancestors all the way to the root, with each “link” in the plot signifying parent-child edge somewhere along the tree sequence) and plot it on a 2D surface (x and y dimensions [0, 1]). The plot is a little chaotic, but hopefully conveys the idea (the “focal node” 0 is highlighted in red). This is essentially the same plot we have in the last figure of our paper.

ancestral_links <- ts_ancestors(ts, 0)

ggplot() +
    geom_sf(data = ancestral_links, size = 0.5, aes(alpha = parent_time)) +
    geom_sf(data = sf::st_set_geometry(ancestral_links, "parent_location"), aes(color = parent_time)) +
    geom_sf(data = data[data$node_id == 0, ], size = 3, color = "red")

Conclusion

In this vignette we gave a brief overview of using slendr’s R-tskit interface for loading, processing, and analyzing “pure” non-slendr tree sequences produced by msprime and SLiM scripts.

Although we have only touched upon the most basic features of its R-tskit interface for standard tree sequences, it is important to note that as far as slendr is concerned, it does not matter how a tree sequence was produced, as long as it conforms to the tskit specification. This means that regardless of the source of your tree sequence data, you should be able to use slendr’s tskit functionality to run your analyses.