Functions

birth_death_shift(model, data)

Calculates average branch rates under the birth-death-shift model with a finite state space.

Example:

using Pesto

phy = readtree(Pesto.path("bears.tre")) 
ρ = 1.0  
data = make_SSEdata(phy, "", ρ; include_traits = false) 
λ = [0.1, 0.2] 
μ = [0.05, 0.15] 

η = 0.05 
model = SSEconstant(λ, μ, η)

res = birth_death_shift(model, data)

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)}\]

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:

ρ = 1.0
λ = 1.0
μ = 0.5
t = 0.1

ψ(t, λ, μ, ρ)

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"))
ρ = 1.0
data = make_SSEdata2(phy, ρ)

lp(λ, μ, data)

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"))
ρ = 0.67
data = make_SSEdata2(phy, ρ)

λ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")) 
ρ = 0.635
data = make_SSEdata2(phy, ρ)

λ = [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 = SSEconstant(λ, μ, ηml)

11

make_descendants(data)

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

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

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"))
ρ = 0.635
data = make_SSEdata(phy, ρ)

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"))
ρ = 0.635
primates = SSEdata(phy, ρ)

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

rates = birth_death_shift(model, primates)

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

make_SSEdata(phy, ρ)

Example:

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

make_descendants_nodes(data)

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

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

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)

postorder_nosave(model::SSEconstant, data::SSEdata, E, alg = OrdinaryDiffEq.Tsit5())

TBW