# This file contains the core functions to calculate the likelihood of the observed gene expression
# given a set of parameters for the EVE model
# Author: Rori Rohlfs, Lars Gronvold, John Mendoza
# Date 2/25/19
#' Calculate expected species mean from phylogeny under the OU model
#'
#' Calculate expected species mean, and evolutionary variance for each node in the tree
#' given theta, sigma^2 and alpha parameters for each edge in the tree.
#'
#' @inheritParams EVEmodel
#'
#' @return Vector of expected species means for each node in the tree
#' @export
#' @importFrom stats cophenetic
#' @import ape
expectedMeanOU <- function(tree, thetas, alphas, rootE)
{
# Pre-allocate vector
expected.mean <- numeric( Ntip(tree) + Nnode(tree) )
# Get the order of the edges to traverse from root to tips
edgeOrder <- rev( postorder(tree) )
# Get index of root
root.index <- tree$edge[edgeOrder[1], 1]
# Set root values
expected.mean[root.index] = rootE
for(i in edgeOrder)
{
parent.index <- tree$edge[i, 1]
child.index <- tree$edge[i, 2]
delta.T <- tree$edge.length[i]
expected.mean[child.index] <- expected.mean[parent.index] * exp(-alphas[i] * delta.T) +
thetas[i] * (1 - exp(-alphas[i] * delta.T))
}
return( expected.mean )
}
#' Calculate evolutionary variance from phylogeny under the OU model
#'
#' Calculate expected species mean, and evolutionary variance for each node in the tree
#' given theta, sigma^2 and alpha parameters for each edge in the tree.
#'
#' @inheritParams EVEmodel
#'
#' @return Vector of variances for each node
#' @export
#' @import ape
evolVarOU <- function(tree, alphas, sigma2s, rootVar)
{
# Pre-allocate vector
evol.variance <- numeric( ape::Ntip(tree) + ape::Nnode(tree) )
# Get the order of the edges to traverse from root to tips
edgeOrder <- rev(ape::postorder(tree))
# Get index of root
root.index <- tree$edge[edgeOrder[1], 1]
# Set root values
evol.variance[root.index] = rootVar
for(i in edgeOrder)
{
parent.index <- tree$edge[i, 1]
child.index <- tree$edge[i, 2]
delta.T <- tree$edge.length[i]
evol.variance[child.index] <- evol.variance[parent.index] * exp(-2 * alphas[i] * delta.T) +
sigma2s[i] / (2 * alphas[i]) * (1 - exp(-2 * alphas[i] * delta.T))
}
return( evol.variance )
}
#' Covariance from phylogeny under the OU model
#'
#' Calculates the OU model covariance matrix given alpha parameters for each edge
#' and the evolutionary variance at each node of a tree.
#'
#' @inheritParams EVEmodel
#' @param evol.variance Vector of evolutionary variance for each edge in the tree (as generated by \code{\link{evolVarOU}})
#'
#' @return Covariance matrix
#' @export
#'
#' @import ape
covarianceOU <- function(tree, alphas, evol.variance)
{
# Copy the phylogeny and multiply edge lengths with alpha
attenuationTree <- tree
attenuationTree$edge.length <- attenuationTree$edge.length * alphas
# calculate the attenuation matrix using the cophenetic distance function
attenuation.Matrix <- stats::cophenetic(attenuationTree)
# get matrix of variances of the most recent common ancestors
variance.MRCA <- apply(ape::mrca(tree), 1:2, function(i) evol.variance[i])
# Return the covariance matrix
return (variance.MRCA * exp(-attenuation.Matrix))
}
# rootVar = beta * sigma2 / (2 * alpha)
# rootE = theta
#' Generate multivariate normal distribution parameters for an EVE model
#'
#' @param tree Species phylogeny (phylo object)
#' @param thetas Vector of theta parameters for each edge in the tree
#' @param alphas Vector of alpha parameters for each edge in the tree
#' @param sigma2s Vector of sigma2 parameters for each edge in the tree
#' @param beta The beta parameter
#' @param rootVar The evolutionary variance at the root node
#' @param rootE The expected species mean at the root node
#' @param index.expand A vector mapping individual samples, i.e. columns in the
#' gene expression matrix, to the species in the phylogeny
#'
#' @return List with covariance matrix (sigma) and expected species means (mean)
#' @export
EVEmodel <- function( tree, thetas, alphas, sigma2s, beta, rootVar, rootE, index.expand ){
E <- expectedMeanOU(tree, thetas, alphas, rootE)
V <- evolVarOU(tree, alphas, sigma2s, rootVar)
sigmaOU <- covarianceOU(tree, alphas, evol.variance = V)
# Expand covariance matrix
sigma <- sigmaOU[index.expand, index.expand]
# Add within-species variance
diag(sigma) <- diag(sigma)*(1+beta)
# Expand expected values
mean <- E[index.expand]
# Return expanded sigma and mean (expected values)
return( list(sigma = sigma, mean = mean) )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.