estimate_ziber: Parameter estimation of zero-inflated bernoulli model

In scone: Single Cell Overview of Normalized Expression data

Description

This function implements an expectation-maximization algorithm for a zero-inflated bernoulli model of transcript detection, modeling gene expression state (off of on) as a bernoulli draw on a gene-specific expression rate (Z in 0,1). Detection conditioned on expression is a logistic function of gene-level features. The bernoulli model is modeled numerically by a logistic model with an intercept.

Usage

estimate_ziber(x, fp_tresh = 0, gfeatM = NULL, bulk_model = FALSE,
  pos_controls = NULL, em_tol = 0.01, maxiter = 100,
  verbose = FALSE)

Arguments

x	matrix. An expression data matrix (genes in rows, cells in columns)
fp_tresh	numeric. Threshold for calling a positive detection (D = 1). Default 0.
gfeatM	matrix. Numeric gene level determinants of drop-out (genes in rows, features in columns)
bulk_model	logical. Use median log-expression of gene in detected fraction as sole gene-level feature. Default FALSE. Ignored if gfeatM is specified.
pos_controls	logical. TRUE for all genes that are known to be expressed in all cells.
em_tol	numeric. Convergence treshold on log-likelihood.
maxiter	numeric. The maximum number of iterations. Default 100.
verbose	logical. Whether or not to print the value of the likelihood at each iteration.

Value

a list with the following elements:

• W coefficients of sample-specific logistic drop-out model

• Alpha intercept and gene-level parameter matrix

• X intercept

• Beta coefficient of gene-specific logistic expression model

• fnr_character the probability, per gene, of P(D=0|E=1)

• p_nodrop 1 - the probability P(drop|Y), useful as weights in weighted PCA

• expected_state the expected value E[Z] (1 = "on")

• loglik the log-likelihood

• convergence 0 if the algorithm converged and 1 if maxiter was reached

Examples

mat <- matrix(rpois(1000, lambda = 3), ncol=10)
mat = mat * matrix(1-rbinom(1000, size = 1, prob = .01), ncol=10)
ziber_out = suppressWarnings(estimate_ziber(mat,
   bulk_model = TRUE,
   pos_controls = 1:10))

scone documentation built on Nov. 8, 2020, 5:20 p.m.