# R/calc_posterior.R In SAVER: Single-Cell RNA-Seq Gene Expression Recovery

#### Documented in calc.post

#' Calculates SAVER posterior
#'
#' Given prediction and prior variance, calculates the Gamma posterior
#' distribution parameters for a single gene.
#'
#' Let \eqn{\alpha} be the shape parameter and \eqn{\beta} be the rate
#' parameter of the prior Gamma distribution. Then, the posterior Gamma
#' distribution will be
#' \deqn{Gamma(y + \alpha, sf + \beta),}
#' where y is the observed gene count and sf is the size factor.
#'
#' @param y A vector of observed gene counts.
#'
#' @param mu A vector of prior means.
#'
#' @param sf Vector of normalized size factors.
#'
#' @param scale.sf Mean of the original size factors.
#'
#' @return A list with the following components
#' \item{\code{estimate}}{Recovered (normalized) expression}
#' \item{\code{se}}{Standard error of expression estimate}
#'
#' @importFrom stats qgamma
#' @export
calc.post <- function(y, mu, sf, scale.sf) {
n <- length(y)
if (length(mu) == 1) {
mu <- rep(mu, n)
}
if (length(sf) == 1) {
sf <- rep(sf, n)
}
if (sum(y) == 0) {
return(list(estimate = rep(0, n), se = rep(0, n)))
}
if (var(mu) == 0) {
prior.beta <- rep(calc.b(y, mu, sf)[1], n)
prior.alpha <- mu*prior.beta
} else{
a <- tryCatch(calc.a(y, mu, sf), error = function(cond) {
return(c(0, Inf))
})
b <- tryCatch(calc.b(y, mu, sf), error = function(cond) {
return(c(0, Inf))
})
k <- tryCatch(calc.k(y, mu, sf), error = function(cond) {
return(c(0, Inf))
})
var.method <- which.min(c(a[2], b[2], k[2]))
if (var.method == 1) {
prior.alpha <- a[1]
prior.beta <- a[1]/mu
} else if (var.method == 2) {
prior.beta <- b[1]
prior.alpha <- mu*b[1]
} else {
prior.alpha <- mu^2/k[1]
prior.beta <- mu/k[1]
}
}
post.alpha <- prior.alpha + y
post.beta <- prior.beta + sf
lambda.hat <- post.alpha/post.beta
se <- sqrt(post.alpha/post.beta^2)
return(list(estimate = unname(ceiling(lambda.hat*1000*scale.sf)/1000),
se = unname(ceiling(se*1000*scale.sf)/1000)))
}


## Try the SAVER package in your browser

Any scripts or data that you put into this service are public.

SAVER documentation built on Nov. 14, 2019, 1:06 a.m.