Generating trees
TreeTools has three ways to generate trees:
- using a birth-death process (i.e. forward simulation);
- using a coalescent process (i.e. backward simulation);
- using a collection of pre-defined simple shapes.
All three methods are accessible through the TreeTools.Generate submodule. For now, this is not exported: most functions below need to be prefaced by an explicit TreeTools.Generate.
Birth-death process
This is rather straightforward:
julia> begin n = 25 # number of lineages to reach for completion b = 1 # birth rate d = 0.1 # death rate end;julia> tree, completed = TreeTools.Generate.birth_death(n, b, d)(Tree{TreeTools.EmptyData}: 49 nodes, 25 leaves, true)julia> length(leaves(tree)) == n || !completedtrue
The second output argument completed is a boolean indicating whether the process reached the target n lineages. If the death rate is non-zero, it is possible that all lineages die before completion.
The implemented process starts with one lineage (the root). Then, if there are $n$ lineages in the tree:
- with rate $b\cdot n$, trigger a birth event: pick a lineage at random and split it in two;
- with rate $d\cdot n$, trigger a death event: stop the lineage and place a leaf at its end;
- if
nreaches the target value for completion or if all lineages are dead, stop the process by placing a leaf at the end of each lineage.
Coalescent process
There are two coalescent models implemented.
- The classical Kingman's process, with a rate of coalescence $\nu = n(n-1)/N$ for $n$ lineages, where $N$ is a parameter (the population size).
- The Yule process, with a rate of coalescence $\nu = b\cdot n$ for $n$ lineages, where $b$ can be interpreted as a birth rate. This is equivalent to a birth only forward process.
Generating a coalescent tree is a two-step process: first construct a Coalescent object, then give it to the function genealogy.
julia> kingman = TreeTools.Generate.KingmanCoalescent(n=25, N=10_000)TreeTools.Generate.KingmanCoalescent(25, 10000.0)julia> yule = TreeTools.Generate.YuleCoalescent(n=5, b=1)TreeTools.Generate.YuleCoalescent(5, 1.0)julia> tree = TreeTools.Generate.genealogy(kingman)Tree{TreeTools.EmptyData}: 49 nodes, 25 leavesjulia> distance(root(tree), first(leaves(tree))) # ~2N on average10279.927265708873julia> tree = TreeTools.Generate.genealogy(yule)____________________________________________ 2 ______________________________| | |____________________________________________ 4 _| | __________________________________________________________ 1 |________________| | ____ 3 |_____________________________________________________| |____ 5julia> distance(root(tree), first(leaves(tree))) # ~b log(n) on average1.3480284800392872
It is possible to create custom coalescent. Below, we create a backward process with multiple mergers, with the following properties:
- the rate of merging is $\rho$, independent of the number of lineages $n$;
- when a merge occurs, an average fraction $\beta$ of the $n$ lineages is involved (sampled using a binomial).
First, we define the corresponding coalescent type. Importantly, the field n must be present in the struct, and is interpreted as the number of remaining lineages.
julia> @kwdef mutable struct MultipleMerger <: TreeTools.Generate.Coalescent n::Int # number of lineages β::Float64 # parameters ρ::Float64 endMain.MultipleMerger
Now, we need to define the choose_event function that samples a merge event from this process. The return value should be a tuple k, t with k being the number of lineages involved in the merge and t the time to the event.
julia> using Distributionsjulia> import TreeTools.Generate.choose_eventjulia> function choose_event(C::MultipleMerger) C.n <= 1 && throw(ArgumentError("Cannot choose coalescence event for $(C.n) lineage.")) t = 0. merger = 0 coin = Binomial(C.n, C.β) # this is from Distributions, while merger < 2 t += rand(Exponential(1/C.ρ)) merger = rand(coin) end return merger, t endchoose_event (generic function with 5 methods)
Now we can sample from this
julia> mm = MultipleMerger(n=15, β=.25, ρ=1)Main.MultipleMerger(15, 0.25, 1.0)julia> tree = TreeTools.Generate.genealogy(mm)_________________________________________________________________________ 2 | | ______________ 6 | | | _________________________________|______________ 7 _| | | | | |______________ 14 | | | | _______________________ 8 | | | |________________________| | __________________ 5 | | | | | |__________________ 13 | |____| | | | ____________ 1 | | | | |________________________| |_____|____________ 4 | | | |____________ 9 | | _____________________ 3 | | | |_____________________ 10 | | |_|_____________________ 11 | |_____________________ 12 | |_____________________ 15julia> using StatsBasejulia> map(n -> length(children(n)), internals(tree)) |> countmap # if we're not unlucky, there should be some multiple mergers in here!Dict{Int64, Int64} with 3 entries: 5 => 1 2 => 2 3 => 4
Basic shapes
There are three at the moment, with self explanatory names:
julia> star = TreeTools.Generate.star_tree(4, 1.)____________________________________________________________________________ 1 | |____________________________________________________________________________ 2 _| |____________________________________________________________________________ 3 | |____________________________________________________________________________ 4julia> ladder = TreeTools.Generate.ladder_tree(4, 1.)____________________________________________________________________________ 4 _| | ___________________________________________________ 3 |________________________| | _________________________ 2 |_________________________| |_________________________ 1julia> balanced = TreeTools.Generate.balanced_binary_tree(4, 1.)______________________________________ 1 _____________________________________| | |______________________________________ 2 _| | ______________________________________ 3 |_____________________________________| |______________________________________ 4
The docstrings give a bit more detail.