R/EVEcore.R

In Jmendo12/evemodel: Expression Variance Evolution Model

# 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) )
}