Simulates the evolution of a quantitative trait.
Functions to simulate the evolution of a quantitative trait along a
phylogenetic tree inputted as an object of class ‘phylo’
(package ape) or graph-class object.
EvolveOptimMarkovTree(tp, tw, anc, p=1, root=tp$edge[1,1]) TraitOUsimTree(tp, a, sigma, opt, p=1, root=tp$edge[1,1]) OUvar(d, a=0, theta=1, sigma=1) PEMvar(d, a=0, psi=1) TraitVarGraphSim(x, variance, distance="distance", p=1, ...)
tp |
A rooted phylogenetic tree of class ‘phylo’ (see package ape), |
tw |
Transition matrix giving the probability that the optimum trait value changes from one state to another at vertices. All rows must sum to 1. |
anc |
Ancestral state of a trait (at root). |
p |
Number of variates to generate. |
root |
Root node of the tree. |
d |
Phylogenetic distances (edge lengths). |
a |
|
theta |
Adaptive evolution rate, i.e. mean trait shift by natural selection. |
sigma |
Neutral evolution rate, i.e. mean trait shift by drift. |
psi |
Mean evolution rate. |
opt |
An index vector of optima at the nodes. |
x |
A |
variance |
Variance function ( |
distance |
The name of the member of ‘x$edge’ where edge lengths can be found. |
... |
Additional parameters for the specified variance function. |
Function EvolveOptimMarkovTree allows one to simulate the changes of
optimum trait values as a Markov process. The index whereby the process
starts, at the tree root, is set by parameter anc; this is the
ancestral character state. From the root onwards to the tips, the
optimum is given the opportunity to change following a multinomial
random draw with transition probabilities given by the rows of matrix
tw. The integers thus obtained can be used as indices of a
vector featuring the actual optimum trait values corresponding to the
simulated selection regimes. The resulting optimum trait values at the
nodes are used by TraitOUsimTree as its parameters
opt to simulate trait values at nodes and tips. Function
TraitVarGraphSim uses a graph variance function (either
OUvar or PEMvar) to reconstruct a covariance matrix that
is used to generate covariates drawn from a multi-normal distribution.
Functions EvolveOptimMarkovTree and
TraitOUsimTree return a matrix whose rows represent the
vertices (nodes and tips) of the phylogenetic tree and whose columns
stand for the n different trials the function was asked to
perform. For EvolveQTraitTree, the elements of the matrix are
integers, representing the selection regimes prevailing at the nodes
and tips, whereas for TraitOUsimTree, the elements are
simulated quantitative trait values at the nodes and tips. These
functions are implemented in C language and therefore run swiftly even
for large (10000+ species) trees.
Function TraitVarGraphSim returns p phylogenetic
signals and is implemented using a rotation of a matrix of standard
normal random (mean=0, variance=1) deviates. The rotation matrix is
itself obtained by Choleski factorization of the trait covariance
matrix expected for a given set of trees, variance function, and
variance function parameters.
Guillaume Guénard, Département de sciences biologiques Université de Montréal, Montréal, QC, Canada.
Butler, M. A. & King, A. A. 2004. Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am. Nat. 164: 683-695.
Guénard, G., Legendre, P., and Peres-Neto, P. 2013. Phylogenetic eigenvector maps (PEM): a framework to model and predict species traits. Meth. Ecol. Evol. In press.
opt <- c(-2,0,2) # Three trait optima: -2, 0, and 2
### Transition probabilities:
transit <- matrix(c(0.7,0.2,0.2,0.2,0.7,0.1,0.1,0.1,0.7),
length(opt),length(opt),dimnames=list(from=opt,to=opt))
# In this example, the trait has a probability of 0.7 to stay at a given
# optimum, a probability of 0.2 for the optimum to change from -2 to 0,
# from 0 to -2, and from 2 to -2, and a probability of 0.1 for the
# optimum to change from -2 to 2, from 0 to 2, and from 2 to 0.
nsp <- 25 # A random tree for 25 species.
tree2 <- rtree(nsp,tip.label=paste("Species",1:nsp,sep=""))
tree2$node.label=paste("N",1:tree2$Nnode,sep="") # Node labels.
## Simulate 10 trials of optimum change.
reg <- EvolveOptimMarkovTree(tp=tree2,tw=transit,p=10,anc=2)
y1 <- TraitOUsimTree(tp=tree2,a=0,sigma=1,
opt=opt[reg[,1]],p=10) # Neutral
y2 <- TraitOUsimTree(tp=tree2,a=1,sigma=1,
opt=opt[reg[,1]],p=10) # Few selection.
y3 <- TraitOUsimTree(tp=tree2,a=10,sigma=1,
opt=opt[reg[,1]],p=10) # Strong selection.
### Display optimum change with colours.
displayOUprocess <- function(tp,trait,regime,mvalue) {
layout(matrix(1:2,1,2))
n <- length(tp$tip.label)
ape::plot.phylo(tp,show.tip.label=TRUE,show.node.label=TRUE,root.edge=FALSE,
direction="rightwards",adj=0,
edge.color=rainbow(length(trait))[regime[tp$edge[,2]]])
plot(y=1:n,x=mvalue[1:n],type="b",xlim=c(-5,5),ylab="",xlab="Trait value",yaxt="n",
bg=rainbow(length(trait))[regime[1:n]],pch=21)
text(trait[regime[1:n]],y=1:n,x=5,col=rainbow(length(trait))[regime[1:n]])
abline(v=0)
}
displayOUprocess(tree2,opt,reg[,1],y1[,1]) # Trait evolve neutrally,
displayOUprocess(tree2,opt,reg[,1],y2[,1]) # under weak selection,
displayOUprocess(tree2,opt,reg[,1],y3[,1]) # under strong selection.
#
x <- Phylo2DirectedGraph(tree2)
y4 <- TraitVarGraphSim(x, variance = MPSEM::OUvar, p=10, a=5)
#
DisplayTreeEvol <- function(tp,mvalue) {
layout(matrix(1:2,1,2))
n <- length(tp$tip.label)
ape::plot.phylo(tp,show.tip.label = TRUE, show.node.label = TRUE, root.edge = FALSE,
direction = "rightwards", adj = 0)
plot(y=1:n,x=mvalue[1:n],type="b",xlim=c(-5,5),ylab="",xlab="Trait value",yaxt="n",pch=21)
abline(v=0)
}
## Recursively displays the simulated traits.
for(i in 1:10) {
DisplayTreeEvol(tree2,y4[i,])
if(is.null(locator(1)))
break # Stops recursive display on a mouse right-click.
}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.