In coleoguy/chromevolR: Analysis of chromosome number evolution and binary traits

The model for ChromPlus

Our model builds on existing dynamic phylogenetic models of chromosome number evolution (mayrose2010, glick2014), and, like the model of zenil 2017, it incorporates a binary trait that can affect rates of chromosome number change. Additionally, our model can allow for the binary trait to affect rates of speciation and extinction. Accounting for the possibility of state dependent diversification is a vital component because of the history of speciation models that invoke chromosomal evolution (bickham1979, grant1981, templeton1981, baker1986). Our model implementation is available in an R package called chromePlus . It is built on the diversitree framework (fitzjohn2012), which provides additional amenities such as allowing for uncertainty in the binary character state, sampling incompleteness, and use of existing MCMC (Markov Chain Monte Carlo) or likelihood maximization functions to estimate model parameters.

A graphical schematic of our model is shown in Figure 1, depicting the ten possible kinds of state transitions. While in binary state 1 a lineage with $i$ chromosomes may first simply stay in its current state. Otherwise, it may either increase chromosome number to $i+1$ (ascending dysploidy, $\lambda_1$), reduce chromosome number to $i-1$ (descending dysploidy, $\delta_1$), or increase chromosome number to $1.5i$ or $2i$ (demi-polyploidy, $\mu_1$ and polyploidy, $\rho_1$ respectively). Likewise, a lineage in binary state 2 with $i$ chromosomes may remain in its current state, increase chromosome number to $i+1$ (ascending dysploidy, $\lambda_2$), reduce chromosome number to $i-1$ (descending dysploidy, $\delta_2$), or increase chromosome number to $1.5i$ or $2i$ (demi-polyploidy, $\mu_2$ and polyploidy, $\rho_2$ respectively). Additionally, a lineage that is currently in one of the binary states can transition to the other state (transition from state 1 to state 2, $q_{12}$ or from state 2 state 1, $q_{21}$); these transitions occur with no change in the chromosome number. Because chromosome number is limited to whole numbers, demi-polyploidy as described above would only be possible for even values of $i$. As in previous work, we assume that demi-polyploidy when $i$ is odd leads to either of the closest whole numbers with equal probability. The state-dependent version of our model has four additional parameters: speciation and extinction in binary state 1, and speciation and extinction in binary state 2.

knitr::include_graphics("model.png", dpi=400)

Figure 1 ChromePlus Model

Simple model chromosome increase, decrease, and whole genome duplication only:

# load packages
library(chromePlus)
library(diversitree)

# read tree
tree <- read.tree("../inst/tree.new")

# read data
dat1 <- read.csv("../inst/tip.data.certain.csv", as.is = T)

# convert data appropriate for analysis
dat.mat <- datatoMatrix(x     = dat1,
                        range = c(18,80),
                        hyper = F)   # this indicates no binary trait

# make the basic likelihood function
lik <- make.mkn(tree = tree,
                states = dat.mat,
                k = ncol(dat.mat),
                strict= F,
                control=list(method="ode"))

# constrain the likelihood function to a biologically realistic design
con.lik <- constrainMkn(data = dat.mat,
                        lik = lik,
                        polyploidy = F,
                        hyper = F,
                        constrain = list(drop.demi = T,
                                         drop.poly= F))

# lets make sure we have the parameters we expect
argnames(con.lik)

Model with a binary state that impacts chromosome evolution

# load packages
library(chromePlus)
library(diversitree)

# read tree
tree <- read.tree("../inst/tree.new")

# read data
dat1 <- read.csv("../inst/tip.data.certain.csv", as.is = T)

# convert data appropriate for analysis
dat.mat <- datatoMatrix(x     = dat1,
                        range = c(18,80),
                        hyper = T)   # this indicates no binary trait

# make the basic likelihood function
lik <- make.mkn(tree = tree,
                states = dat.mat,
                k = ncol(dat.mat),
                strict= F,
                control=list(method="ode"))

# constrain the likelihood function to a biologically realistic design
con.lik <- constrainMkn(data = dat.mat,
                        lik = lik,
                        polyploidy = F,
                        hyper = T,
                        constrain = list(drop.demi = T,
                                         drop.poly= F))

# lets make sure we have the parameters we expect
argnames(con.lik)

More complex model with uncertainty in tip states, different rates of chromosomes evolution, speciation, and extinction rates associated with the binary character.

# load packages
library(chromePlus)
library(diversitree)

# read tree
tree <- read.tree("../inst/tree.new")

# read data
dat2 <- read.csv("../inst/tip.data.uncertain.csv", as.is = T)

# prepare data
dat.mat2 <- datatoMatrix(x     = dat2,
                         range = c(18,80),
                         hyper = T)

# make the basic likelihood function
lik2 <- make.musse(tree = tree,
                   states = dat.mat2,
                   k = ncol(dat.mat2),
                   strict= F,
                   control=list(method="ode"))

# constrain to biological reality
con.lik2 <- constrainMuSSE(data = dat.mat2,
                           lik = lik2,
                           s.lambda = F,
                           s.mu = F,
                           polyploidy = F,
                           hyper = T,
                           constrain = list(drop.demi = T,
                                            drop.poly= F))
argnames(con.lik2)

In all cases once we have a biologically realistic model we can fit it using find.mle or mcmc from the package diversitree. A brief example is shown below in depth discussion of these functions is available in diversitree.

find.mle(con.lik2, x.init = runif(min = 0, max = 1, n = 5))
mcmc(con.lik2, x.init = runif(min = 0, max = 1, n = 5), w = 1, nsteps = 10)

constraining models

The function constrainMkn has a series of arguments that allow for realistic models of chromosome evolution. The default settings for the functions is to fit the most complex model with all possible transitions shown in figure one estimated seperately. The table below details the arguments that can be used to simplify this full model.

Argument |Effect ------------|------------------------------------------ hyper |when set to FALSE the binary state is ignored so only a single rate of gain, loss, demiploidy, and polyploidy are possible polyploidy |if set to TRUE the the hyper state that is evaluated is polyploidy this means that transition from state 1 to state 2 in the binary character also double the chromosome number constrain |is given a list with any of the 5 additional constraint descriptions given below | drop.poly |if set to TRUE then the rate of polyploidy is set to zero drop.demi |if set to TRUE then the rate of demiploidy is set to zero symmetric |if set to TRUE then the rates estimated in the two binary states are forced to be equal nometa |if set to TRUE then the binary state is ignored meta="ARD" |can be set to "ARD" or "SYM" to describe wether q12 and q21 are the same or different

The function constrainMuSSE has a series of arguments that allow for realistic models of chromosome evolution. The default settings for the functions is to fit the most complex model with all possible transitions shown in figure one estimated seperately. The table below details the arguments that can be used to simplify this full model.

Argument |Effect ------------|------------------------------------------ hyper |when set to FALSE the binary state is ignored so only a single rate of gain, loss, demiploidy, and polyploidy are possible polyploidy |if set to TRUE the the hyper state that is evaluated is polyploidy this means that transition from state 1 to state 2 in the binary character also double the chromosome number s.lambda |if set to TRUE then a single lambda (speciation rate) is estimated s.mu |if set to TRUE then a single mu (extinction rate) is estimated constrain |is given a list with any of the 5 additional constraint descriptions given below | drop.poly |if set to TRUE then the rate of polyploidy is set to zero drop.demi |if set to TRUE then the rate of demiploidy is set to zero symmetric |if set to TRUE then the rates estimated in the two binary states are forced to be equal nometa |if set to TRUE then the binary state is ignored meta="ARD" |can be set to "ARD" or "SYM" to describe wether q12 and q21 are the same or different

coleoguy/chromevolR documentation built on July 27, 2023, 12:40 p.m.