Simple analysis
Here is an example of an analysis of branch-specific rates under the birth-death-shift model.
Tree file
First, we load the necessary modules and read in the tree file.
using Pesto
phy = readtree(Pesto.path("primates.tre"))
ρ = 0.635
primates = SSEdata(phy, ρ)Analysis
A simple analysis can be done like so:
model, rates = pesto(primates)To see how this analysis is set up, see the next section (Extended analysis).
Tree plots
If we want to plot the results immediately, we will need to use the Makie dependency. Makie has several backends, but we will use CairoMakie for now, as it is designed to plot high-quality 2-dimensional figures.
using Makie, CairoMakie
treeplot(primates, rates)
If you want to instead plot the number of shifts, and for example use a different color gradient, you can enter the following
cmap = Makie.cgrad([:black, :red])
f = treeplot(primates, rates, "nshift"; cmap = cmap)
The figure can be saved using the following code
Makie.save("path/to/figure.pdf", f)See section Plot with ggtree for more instructions on how to make more customized and publication-quality figures.
Save to file
A summary of the results of the analysis is stored in the rates object, which is a data frame. We can print the first five rows to get an idea of it:
rates[1:5,:]| Row | mean_lambda | mean_mu | mean_netdiv | mean_relext | delta_lambda | delta_mu | delta_netdiv | delta_relext | node | edge | nshift | shift_bf | shift_bf_log |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Int64 | Int64 | Float64 | Float64 | Float64 | |
| 1 | 0.219675 | 0.137233 | 0.0824416 | 0.614623 | -0.000269656 | -0.000376738 | 0.000107082 | -0.000749151 | 235 | 1 | 0.00656303 | 0.147156 | -0.832222 |
| 2 | 0.219822 | 0.137151 | 0.0826702 | 0.613974 | 0.000750585 | 9.25196e-5 | 0.000658066 | -0.00152837 | 236 | 2 | 0.0200711 | 0.187898 | -0.726078 |
| 3 | 0.221556 | 0.137215 | 0.0843406 | 0.610394 | 0.00324677 | 0.000146542 | 0.00310023 | -0.00601478 | 237 | 3 | 0.0298147 | 0.722026 | -0.141447 |
| 4 | 0.224178 | 0.137536 | 0.0866416 | 0.606666 | 0.00152271 | 0.000592671 | 0.000930041 | -8.00436e-5 | 238 | 4 | 0.0108046 | 0.66752 | -0.175535 |
| 5 | 0.225484 | 0.138266 | 0.087218 | 0.607953 | 0.000992493 | 0.000933741 | 5.87518e-5 | 0.003101 | 1 | 5 | 0.00982901 | 0.813542 | -0.0896202 |
Each row in the data frame represents a branch in the tree, indexed by the edge index edge. Alternatively, you can index the rows using the node index. It is possible to save the data-frame as is to a file, using using CSV; CSV.write("rates.csv", rates). We can also represent the output as an extended newick string:
newick(primates, rates)"(((((Galago_matschiei[&mean_lambda=0.22548398910066372,mean_mu=0.13826595822387158,mean_netdiv=0.08721803087679214,mean_relext=0.6079525096588052,delta_lambda=0.0009924929387501447,delta_mu=0.0009337411730851175,delta_netdiv=5.875176566508267e-5,delta_relext=0.00310100" ⋯ 175205 bytes ⋯ "xt=0.614939107379384,delta_lambda=1.702967027114588e-5,delta_mu=1.957572879640601e-5,delta_netdiv=-2.5460585252184975e-6,delta_relext=2.4161803716671315e-5,edge=94,nshift=0.0010637616944320911,shift_bf=0.1284896504067899,shift_bf_log=-0.8911318525057513]:2.491564):0.0;"If we want to save the newick string to a file, we can use the writenewick function
writenewick("primates_analysis.tre", primates, rates)This tree file can be loaded in other programs such as R and can be plotted using standard packages like ape and ggtree (see section Plot with ggtree).