1 Simulating trees
In this section, we’ll simulate ultrametric and mutation-based trees. Ultrametric just means that the tips are all the same distance from the root. In our specific case, the ultrametric trees will have edge lengths in units of time, so ultrametric means that the tips are sampled at the same time. Mutation-based trees are those with edge lengths in units of mutations. Mutation-based trees will typically not be ultrametric, because there will be some fluctuations in the number of mutations acquired.
We’ll be generating trees of class phylo, which is a
fairly straightforward encoding of phylogenetic trees in R. The ape
documentation describing the class phylo can be found here.
1.1 Ultrametric trees
Let’s start with ultrametric trees. First, let’s simulate a single
tree and have a look at the class phylo. We’ll use the
function simUltra() from our package to do this, which
simulates an ultrametric tree. We’ll set a birth rate a, a
death rate b, an age for the tree cloneAge,
and the number of sampled tips n. Note that we use
cloneAge instead of \(T\)
for the time because T in R means TRUE.
# Generate the tree
tree <- simUltra(a = 1, b = 0, cloneAge = 20, n = 100, addStem = TRUE)
# We see that the tree is of class phylo
class(tree)
#> [1] "phylo"
# Preview the tree
print(tree)
#>
#> Phylogenetic tree with 100 tips and 100 internal nodes.
#>
#> Tip labels:
#> t8, t82, t63, t21, t19, t72, ...
#>
#> Rooted; includes branch lengths.Again, the best place to get more info about the class
phylo is here.
Our functions for simulating trees will include metadata. Let’s take a
look at the metadata included in the tree.
We see that the columns are:
r,a,b: The growth rate params.ris the net growth rate,ais birth rate andbis death rate, sor=a-bcloneAge: The time that passes, in the same units asr,a, andb.n: The number of samples from the birth-death process, producing a tree withntips.runtime_seconds: The elapsed time to generate the tree.addStem: Tells us whether the tree has a root edge preceding the first split. We’ll show the “stem” in plots below.
Now let’s plot the tree.
plot.phylo(tree, direction = "downwards", show.tip.label = FALSE)
axisPhylo(side = 2, backward = FALSE, las = 1)
title(main = "Simulated tree", ylab = "Time")
In the tree we plotted above, we see that the “stem” is the edge extending from time 0 to the first split. We call these splits coalescence events and they’ll be very important for our methods of growth rate estimation.
Also, we see that the tree is ultrametric, meaning that they all are sampled at the same time, in this case 20 units of time after the birth-death process began. As a default, we’ll talk about units of time in years, but that is arbitrary. Just make sure the units of the growth rate (per year) are the same as the time units (years).
Let’s apply our growth rate estimates to this tree, seeing if our
estimates are close to the value of r that we used to
generate the tree. We have two functions for estimating the growth rate
of an ultrametric tree, and we’ll compare their performance later
on:
internalLengths(): This function sums the edge lengths of the internal edges of the tree (excluding the stem). We call this sum \(L_i\). The growth rate is calculated as \(r = L_i / n\) where \(n\) is the number of sampled tips.maxLikelihood(): This function uses a maximum likelihood estimation of the distribution of coalescence times, using the information about when the tree splits. The distribution of coalescence times is approximated by a standard logistic distribution scaled by \(1/r\). Maximizing the likelihood of the observed coalescence times for the \(r\) term gives us our estimate.
Both methods are detailed in our paper.
We know that our actual growth rate is 1, so we can evaluate how our
estimates are doing. One estimate doesn’t really tell us much though.
Let’s try 100. To do this, we just set the nTrees param
equal to 100 in our simUltra() function. You can
parallelize this process if you have the parallel package
installed by simply setting nCores, but we’ll leave this at
1 for now. By printing the elapsed time, we see that with only one core
it takes about a minute.
ptm <- proc.time()
tree.list <- simUltra(a = 1, b = 0, cloneAge = 20, n = 100, nTrees = 100, nCores = 1, addStem = TRUE)
print(proc.time()[["elapsed"]] - ptm[["elapsed"]])
#> [1] 149.59Let’s apply our methods to these trees. We can input either a single
phylo object, or a list of phylo objects to
our growth rate functions. Either way, it’ll return a
data.frame which each row corresponding to an estimate for
each tree.
We gave each method 100 trees to estimate, let’s plot these estimates. Remember that all the input trees had a growth rate of 1, so we want these estimates to be close to 1. We’ll use ggplot’s density plot for this:
# Combine for ggplot formatting
resultsCombined <- rbind(resultsLengths, resultsMaxLike)
# Plot, adding a vertical line at r=1 because that's the true growth rate
ggplot(resultsCombined) +
geom_density(aes(x = estimate, color = method), linewidth = 1.5) +
geom_vline(xintercept = exampleUltraTrees[[1]]$metadata$r) +
theme_bw() +
theme(
axis.text.y = element_blank(), axis.ticks.y = element_blank(),
legend.title = element_blank()
) +
xlab("Net growth rate estimate (r)") +
ylab("Density") +
scale_color_manual(labels = c("Internal lengths", "Max. likelihood"), values = c("black", "#009E73"))
Not bad, but we’ll want to test this more formally over a large range of parameter values. More on that later in this vignette. For now, let’s shift our attention to mutation trees.
1.2 Mutation trees
There are two ways to generate mutation trees:
- We can add mutations to a time-based ultrametric tree
- We can generate a completely new mutation-based tree
Either way, the process of generating mutation trees consists of:
- Generating an ultrametric, time-based tree. This step is only required to generate a new tree
- For each edge, draw the number of mutations from a poisson
distribution with mean equal to the mutation rate
numultiplied by the edge length of the time-based tree.
Let’s take the ultrametric trees we just generated and convert them
to mutation-based trees. For this, we provide the function
ultra2mut():
# Set mutation rate equal to 10 muts/year for all trees
mutTree.list <- ultra2mut(tree.list, nu = 10)Alternatively, we can generate a new set of mutation trees using the
function simMut(). We won’t run this because it will take
another few minutes, and we already have the mutation-based trees
generated above.
# Set params for ultra tree + a mutation rate
mutTree.list2 <- simMut(a = 1, b = 0, cloneAge = 20, n = 100, nTrees = 100, nu = 100, nCores = 1)If we want to estimate the growth rate from a mutation tree, we use
the sharedMuts() function, which works very similarly to
the internalLengths() function. However, instead of
counting the internal edge lengths in units of time, the
sharedMuts() function counts the internal/shared mutations
and uses the mutation rate to scale to time. So, if \(M_i\) is the number of internal or shared
mutations, \(\nu\) (nu) is
the mutation rate, and \(n\) is the
number of samples, the growth rate estimate is \(r = M_i /(\nu n)\).
Let’s apply the sharedMuts() function to the trees that
we converted to mutation-based. Remember that these mutation-based trees
are generated from the ultrametric trees we applied our functions
maxLikelihood() and internalLengths() to.
Let’s plot the estimates all together. First, we’ll combine the
data.frames produced. The sharedMuts() function outputs an
additional column corresponding to the mutation rate nu.
Because of this, using rbind on the full data.frames will give an error.
We only want to look at the estimates for now, so we’ll just combine the
necessary columns from each data.frame.
# Combine the columns with the estimates
colsUse <- c("lowerBound", "estimate", "upperBound", "method")
resultsAll <- rbind(resultsShared[, colsUse], resultsCombined[, colsUse])
# Plot, adding a vertical line at r=1 because that's the true growth rate
ggplot(resultsAll) +
geom_density(aes(x = estimate, color = method), linewidth = 1.5) +
geom_vline(xintercept = exampleUltraTrees[[1]]$metadata$r) +
theme_bw() +
theme(
axis.text.y = element_blank(), axis.ticks.y = element_blank(),
legend.title = element_blank()
) +
xlab("Net growth rate estimate (r)") +
ylab("Density") +
scale_color_manual(labels = c("Internal lengths", "Max. likelihood", "Shared Muts."), values = c(colorPal[1], colorPal[4], colorPal[6]))
It looks like we’re on the right track. From this, it’s hard to tell which is best. Maximum likelihood seems a bit biased on the high side, but the estimates seem to fall in a more narrow range. Let’s do a more quantitative analysis now that we have the hang of it.
