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 average20931.794549823353julia> tree = TreeTools.Generate.genealogy(yule)________________________ 1 __________________________________________________| | | ,_ 2 | |______________________| _| |_ 4 | | _______________________________________________ 3 |___________________________| |_______________________________________________ 5julia> distance(root(tree), first(leaves(tree))) # ~b log(n) on average0.37721699073515935
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)__________________ 1 ______________________________________________________| | | _________ 2 | |________| | |_________ 8 _| | _________________________________ 10 | | | | _________ 6 | | | |_______________________________________| |_________ 9 | ____________| | | |_________ 15 | | | | | | ______ 11 |__________| |__| | |______ 13 | | ______________ 5 | | |_______| ____________ 4 | | | |____________ 7 |_| | __ 3 | | |_________|__ 12 | |__ 14julia> 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: 4 => 1 2 => 7 3 => 2
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.)______________________________________ 3 _____________________________________| | |______________________________________ 4 _| | ______________________________________ 1 |_____________________________________| |______________________________________ 2
The docstrings give a bit more detail.