Functions

birth_death_shift(model, data)

Calculates average branch rates under the birth-death-shift model with a discrete rate categories.

Example:

using Pesto

phy = readtree(Pesto.path("bears.tre")) 
sampling_probability = 1.0  
bears = SSEdata(phy, sampling_probability)

λ = [0.1, 0.2] 
μ = [0.05, 0.12] 
η = 0.05 
model = BDSconstant(λ, μ, η)

rates = birth_death_shift(model, bears)

birth_death_shift(model, data)

Calculates average branch rates under the birth-death-shift model with a discrete rate categories.

Example:

using Pesto

phy = readtree(Pesto.path("bears.tre")) 
sampling_probability = 1.0  
bears = construct_tree(phy, sampling_probability)

λ = [0.1, 0.2] 
μ = [0.05, 0.12] 
η = 0.05 
model = BDSconstant(λ, μ, η)

rates = birth_death_shift(model, bears)

estimate_constant_bdp(data::SSEdata[; xinit = [0.11, 0.09], lower = [0.0001, 0.0001], upper = [20.0, 20.0]])

Estimates the speciation and extinction rate under the reconstructed birth-death process with time-homogeneous rates.

Example:

phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.67
data = make_SSEdata(phy, sampling_probability)

λml, μml = estimate_constant_bdp(data)

optimize_eta(λ, μ, data)

Finds the maximum-likelihood parameter value for η (the transition rate) under the birth-death-shift model with a finite state space, conditional on λ and μ.

Example:

using Pesto

phy = readtree(Pesto.path("primates.tre")) 
sampling_probability = 0.635
data = SSEdata(phy, sampling_probability)

λ = [0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.1, 0.2, 0.3] 
μ = [0.09, 0.09, 0.09, 0.19, 0.19, 0.19, 0.29, 0.29, 0.29] 

ηml = optimize_eta(λ, μ, data)

#model = BDSconstant(λ, μ, ηml)

11

make_descendants(data)

takes node indices as input, and returns edge indices Example:

using Pesto
phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.635
data = make_SSEdata(phy, sampling_probability)

make_descendants(data)

with result

Dict{Int64, Vector{Any}} with 232 entries:
  402 => [330, 331]
  413 => [357, 360]
  425 => [380, 381]
  ⋮   => ⋮

make_ancestors(data)

takes node indices as input, and returns edge indices

Example:

using Pesto
phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.635
data = make_SSEdata(phy, sampling_probability)

make_ancestors(data)

with result

Dict{Int64, Int64} with 464 entries:
  56  => 115
  35  => 71
  425 => 379
  ⋮   => ⋮

lrange(from, to, length)

Similar to range, but with proportional spacing.

Example:

using Pesto

lrange(0.001, 100.0, 6)

with result

6-element Vector{Float64}:
   0.0010000000000000002
   0.010000000000000004
   0.10000000000000002
   1.0000000000000004
  10.000000000000002
 100.00000000000004

allpairwise(λ, μ)

Example:

using Pesto

lambda = [0.2, 0.3]
mu = [0.05, 0.10, 0.15, 0.20] 

λ, μ = allpairwise(lambda, mu)

with result

([0.2, 0.3, 0.2, 0.3, 0.2, 0.3, 0.2, 0.3], [0.05, 0.05, 0.1, 0.1, 0.15, 0.15, 0.2, 0.2])

writenewick(filename, data, rates)

writes a newick file with the rate values as comments

Example:

using Pesto

phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.635
primates = SSEdata(phy, sampling_probability)

λ = [0.1, 0.2, 0.3, 0.4, 0.20]
μ = [0.05, 0.15, 0.05, 0.15, 0.25]
η = 1 / tree_length

model = BDSconstant(λ, μ, η)

rates = birth_death_shift(model, primates)

writenewick("/tmp/newick.tre", primates, rates)

SSEdata(phy, sampling_probability)

Example:

using Pesto
phy = readtree(Pesto.path("primates.tre")) 
sampling_probability = 0.635  
data = SSEdata(phy, sampling_probability) 

make_descendants_nodes(data)

takes node indices as input, and returns node indices Example:

using Pesto
phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.635
data = make_SSEdata(phy, sampling_probability)

make_descendants_nodes(data)

readtree("/path/to/phylo.tre")

reads a NEXUS or Newick file using RCall and the R-package ape

Example:

using Pesto
phy = readtree(Pesto.path("primates.tre"))

display(phy)

lp(λ, μ, data)

From Louca and Pennell 2020 (Nature), eq. S28

\[L = \frac{\rho^{n+1}} {\lambda (1 - E(t_1))^2} \times \prod_{i=1}^n \lambda \times \psi(0, t_i) \\ E(t) = 1 - \frac{\exp(\lambda - \mu)t}{\frac{1}{\rho} + \frac{\lambda}{\lambda -\mu} \Big ( \exp((\lambda - \mu)t) - 1 \Big)} \\ \psi(t) = \frac{e^{t(\lambda - \mu)}}{ [ 1 + \frac{\rho \lambda}{\lambda - \mu}(e^{t(\lambda - \mu)} - 1)]^{2}}\]

Logged:

\[\log(L) = (n+1) \log(\rho) + \log(\psi(t_1)) - \log(\lambda) - 2 \log(1 - E(t_1)) + \sum_{i=1}^n \log(\lambda) + \log(\psi(t_i))\]

Example:

λ = 1.0
μ = 0.5

phy = readtrees(Pesto.path("bears.tre"))
sampling_probability = 1.0
data = make_SSEdata(phy, sampling_probability)

lp(λ, μ, data)

Equation S5 in Morlon et al. 2011 [PNAS]

\[\psi(s, t) = e^{(\lambda - \mu)(t - s)} [ 1 + \frac{\frac{\lambda}{\lambda - \mu}(e^{t(\lambda - \mu)} - e^{s(\lambda-\mu)})}{\frac{1}{\rho} + \frac{\lambda}{\lambda - \mu} \times (e^{s(\lambda-\mu)}-1)}]^{-2}\]

We use this one, simplified where s = 0

\[\psi(t) = \frac{e^{t(\lambda - \mu)}}{ [ 1 + \frac{\rho \lambda}{\lambda - \mu}(e^{t(\lambda - \mu)} - 1)]^{2}}\]

Example:

sampling_probability = 1.0
λ = 1.0
μ = 0.5
t = 0.1

psi(t, λ, μ, sampling_probability)

from Morlon et al. 2011 [PNAS], eq. S4

\[E(t) = 1 - \frac{\exp(t(\lambda - \mu))}{\frac{1}{\rho} + \frac{\lambda}{\lambda -\mu} \Big ( \exp((\lambda - \mu)t) - 1 \Big)}\]

magnitude(model, data)

Computes the overall point estimate of the magnitude of the rate shifts in the tree, for a given model. The equation for the magnitude is

\[\text{mag} = \frac{1}{\sum_j \sum_i \hat{N}_{ij}} \sum_j \sum_i (r_i - r_j) \hat{N}_{ij},\]

where $\hat{N}_{ij}$ is the number of estimated diversification rate shifts (summed across all branches) that depart from the rate category j and arrive in the rate category i (in the direction of old to young).

Example:

phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.635
primates = SSEdata(phy, sampling_probability)

λ = [0.05, 0.15, 0.25]
μ = [0.03, 0.08, 0.05]
η = 0.003

model = BDSconstant(λ, μ, η)

mag = magnitude(model, primates)
println(mag)

postorder_async(model, tree)

This function does the postorder pass, and returns a tuple of (u, sf). The matrix u is of size K x 2, where K are the number of rate categories.

The first column in u represents the extinction probability P(Ψ goes extinct|Z == j, θ), conditional on that the rate category was j at the root, and conditional on the parameters of the model θ (i.e. all of the diversification rate categories λj, μj, ψj, as well as the sampling probability ρ).

The second column Dj represents Dj = P(Ψ|Z == j,θ), i.e. the probability density of observing the reconstructed tree Ψ, conditional on the same as above.

The scaling factor sf is there because we normalize the probability densities Dj such that they sum to one.

Example:

λ = [0.1, 0.2, 0.3]
μ = [0.05, 0.05, 0.05]
η = 0.01

model = BhDhModel(λ, μ, η)

phy = readtree(Pesto.path("bears.tre"))
sampling_probability = 1.0
tree = construct_tree(phy, sampling_probability)

Ds, sf = postorder_async(model, tree)

pesto_twostep(data)

this function performs a two-step analysis, where the estimation of the parameters is split in two steps:

step one

• we set up two variables ($\hat{\lambda}$ and $\hat{\mu}$) representing the mean birth and death rates, and we find an estimate for them by maximizing the likelihood of the parameters under the lineage-homogeneous birth-death process

step two

• we set up a birth-death-shift model where the means of the base distributions for the speciation and extinction rates are $\hat{\lambda}$ and $\hat{\mu}$, and the shift rate (η) is a parameter
• we find the estimate for η by maximizing the likelihood under the birth-death-shift process

branch results

finally, we compute branch-specific results, including i) mean rates, ii) mean no. rate shifts, iii) Bayes factors, given the model found by the two-step procedure

Example:

using Pesto

phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.635
data = SSEdata(phy, sampling_probability)

model, rates = pesto_twostep(data)

pesto_fossil(tree)

similar to pesto(...), however for a fossilized birth-death-shift (FBDS) process.

This function does the following

1. sets up a FBDS model where the speciation rate and fossilization rate is allowed to vary, but extinction rate is assumed to be constant
2. estimates five parameters using maximum likelihood: i) the mean speciation rate, ii) the mean fossilization rate, iii) the (constant) extinction rate), iv) the speciation shift rate, and v) the fossilization shift rate.
3. compute branch-specific results given the model with the ML estimates

Example:

using Pesto

phy = readtree(Pesto.path("osteoglossomorpha.tre"))
sampling_probability = 0.234
tree = construct_tree(phy, sampling_probability)

## replace extant leaves that are at least 0.5 Ma older
## than the present with terminal fossil sampling events
## (with survival probability = 1)
assign_fossils!(tree, 0.5)

model, rates = pesto_fossil(tree)

pesto(data)

performs maximum-likelihood estimation of the base distribution parameter and the shift rate, and calculates branch-specific posterior mean i) rates, ii) no. shifts and iii) Bayes factors

Example:

using Pesto

phy = readtree(Pesto.path("primates.tre"))
sampling_probability = 0.635
data = SSEdata(phy, sampling_probability)

model, rates = pesto(data)