find.lineages {paleobuddy} | R Documentation |
Separate a paleobuddy simulation into monophyletic clades
Description
Separates a sim
object into sim
objects each with a mother
species and its descendants. If argument S
is not used, it returns by
default the list of sim
objects descended from each species with an
NA
parent in the original input (meaning species alive at the
beginning of the simulation). If a vector of numbers is supplied for
S
, the list of sim
objects return will instead be descended
from each species in S
. Returns for each clade a vector with the
original identity of member species as well.
Usage
find.lineages(sim, S = NULL)
Arguments
sim |
A |
S |
A vector of species in |
Value
A list
object with (named) sim
objects corresponding
to the clades descended from species in S
. For each clade, an extra
vector LIN
is included so the user can identify the order of species
in the returned sim
objects with the order of species in the original
simulation.
Author(s)
Bruno do Rosario Petrucci and Matheus Januario.
Examples
###
# first, we run a simple simulation with one starting species
# set seed
set.seed(1)
# run simulation with a minimum of 20 species
sim <- bd.sim(n0 = 3, lambda = 0.1, mu = 0.1, tMax = 10,
nFinal = c(20, Inf))
# get a simulation object with the clade originating from species 2
clades <- find.lineages(sim, S = 2)
# now we can check to make sure the subclade was correctly separated
# change NA to 0 on the clade's TE
clades[[1]]$sim$TE[clades[[1]]$sim$EXTANT] <- 0
# plot the phylogeny
if (requireNamespace("ape", quietly = TRUE)) {
plot <- ape::plot.phylo(
make.phylo(clades[[1]]$sim),
main = "red: extinction events \n blue: speciation events");
ape::axisPhylo()
}
# check speciation times
for (j in 2:length(clades[[1]]$sim$TS)) {
# the subtraction is just to adjust the wt with the plot scale
lines(x = c(
sort(clades[[1]]$sim$TS, decreasing = TRUE)[2] -
clades[[1]]$sim$TS[j],
sort(clades[[1]]$sim$TS, decreasing = TRUE)[2] -
clades[[1]]$sim$TS[j]),
y = c(plot$y.lim[1], plot$y.lim[2]), lwd = 2, col = "blue")
}
# check extinction times:
for (j in 1:length(sim$TE)) {
# the subtraction is just to adjust the wt with the plot scale
lines(x = c(
sort(clades[[1]]$sim$TS, decreasing = TRUE)[2] -
clades[[1]]$sim$TE[j],
sort(clades[[1]]$sim$TS, decreasing = TRUE)[2] -
clades[[1]]$sim$TE[j]),
y = c(plot$y.lim[1], plot$y.lim[2]), lwd = 2, col = "red")
}
###
# now we try a simulation with 3 clades
# set seed
set.seed(4)
# run simulation
sim <- bd.sim(n0 = 3, lambda = 0.1, mu = 0.1, tMax = 10,
nFinal = c(20, Inf))
# get subclades descended from original species
clades <- find.lineages(sim)
# get current par options so we can reset later
oldPar <- par(no.readonly = TRUE)
# set up for plotting side by side
par(mfrow = c(1, length(clades)))
# for each clade
for (i in 1:length(clades)) {
# change NA to 0 on the clade's TE
clades[[i]]$sim$TE[clades[[i]]$sim$EXTANT] <- 0
# if there is only one lineage in the clade, nothing happens
if (length(clades[[i]]$sim$TE) < 2) {
# placeholder plot
plot(NA, xlim = c(-1, 1), ylim = c(-1, 1))
text("simulation with \n just one lineage", x = 0, y = 0.5, cex = 2)
}
# else, plot phylogeny
else {
if (requireNamespace("ape", quietly = TRUE)) {
plot <- ape::plot.phylo(
make.phylo(clades[[i]]$sim),
main = "red: extinction events \n blue: speciation events");
ape::axisPhylo()
}
# check speciation times
for (j in 2:length(clades[[i]]$sim$TS)) {
# the subtraction is just to adjust the wt with the plot scale
lines(x = c(
sort(clades[[i]]$sim$TS, decreasing = TRUE)[2] -
clades[[i]]$sim$TS[j],
sort(clades[[i]]$sim$TS, decreasing = TRUE)[2] -
clades[[i]]$sim$TS[j]),
y = c(plot$y.lim[1], plot$y.lim[2]), lwd = 2, col = "blue")
}
# check extinction times:
for (j in 1:length(sim$TE)) {
# the subtraction is just to adjust the wt with the plot scale
lines(x = c(
sort(clades[[i]]$sim$TS, decreasing = TRUE)[2] -
clades[[i]]$sim$TE[j],
sort(clades[[i]]$sim$TS, decreasing = TRUE)[2] -
clades[[i]]$sim$TE[j]),
y = c(plot$y.lim[1], plot$y.lim[2]), lwd = 2, col = "red")
}
}
}
# reset par
par(oldPar)
###
# we can also have an example with more non-starting species in S
# set seed
set.seed(3)
# run simulation
sim <- bd.sim(n0 = 1, lambda = 0.1, mu = 0.1, tMax = 10,
nFinal = c(10, Inf))
# get current par options so we can reset later
oldPar <- par(no.readonly = TRUE)
# set up for plotting side by side
par(mfrow = c(1, 2))
if (requireNamespace("ape", quietly = TRUE)) {
# first we plot the clade started by 1
ape::plot.phylo(make.phylo(sim), main = "original")
ape::axisPhylo()
# this should look the same
ape::plot.phylo(make.phylo(find.lineages(sim)[[1]]$sim),
main="after find.lineages()")
ape::axisPhylo()
# get sublcades descended from the second and third species
clades <- find.lineages(sim, c(2,3))
# and these should be part of the previous phylogenies
ape::plot.phylo(make.phylo(clades$clade_2$sim),
main = "Daughters of sp 2")
ape::axisPhylo()
ape::plot.phylo(make.phylo(clades$clade_3$sim),
main = "Daughters of sp 3")
ape::axisPhylo()
}
# reset par
par(oldPar)
###
# if there is only one clade and we use the default for
# S, we get back the original simulation object
# set seed
set.seed(1)
# run simulation
sim <- bd.sim(n0 = 1, lambda = 0.1, mu = 0.08, tMax = 10,
nFinal = c(5, Inf))
# get current par options so we can reset later
oldPar <- par(no.readonly = TRUE)
# set up for plotting side by side
par(mfrow = c(1, 2))
# plotting sim and find.lineages(sim) - should be equal
if (requireNamespace("ape", quietly = TRUE)) {
ape::plot.phylo(make.phylo(sim), main="original")
ape::axisPhylo()
ape::plot.phylo(make.phylo(find.lineages(sim)[[1]]$sim),
main="after find.lineages()")
ape::axisPhylo()
}
# reset par
par(oldPar)