SIMPLE: SIMPLE: Imputing zero entries and clustering for scRNASeq...

In JunLiuLab/SIMPLEs: SIMPLEs: single-cell RNA sequencing imputation and cell clustering methods by modeling gene module variation

Description

SIMPLE imputes zeros in the gene expression data using the expression level in similar cells and gene-gene correlation. Zero entries in the observed expression matrix come from molecule loss during the experiment ('dropout') or too low expression to be measured. We used Monte Carlo EM algorithm to sample the imputed values and obtain MLEs of all the parameters.

Usage

SIMPLE(
  dat,
  K0 = 10,
  M0 = 1,
  iter = 10,
  est_lam = 1,
  impt_it = 5,
  penl = 1,
  init_imp = NULL,
  sigma0 = 100,
  pi_alpha = 1,
  beta = NULL,
  verbose = F,
  max_lambda = T,
  lambda = NULL,
  sigma = NULL,
  mu = NULL,
  est_z = 1,
  clus = NULL,
  p_min = 0.4,
  cutoff = 0.1,
  K = 20,
  min_gene = 2000,
  num_mc = 3,
  fix_num = F,
  mcmc = 50,
  burnin = 2
)

Arguments

dat	scRNASeq data matrix. Each row is a gene, each column is a cell.
K0	Number of latent gene modules. Default is 10. See details.
M0	Number of clusters. Default is 1. See details.
iter	The number of EM iterations using full data set. See details.
est_lam	The iteration starts to estimate lambda.
impt_it	The iteration starts to sample new imputed values in initial phase. See details.
penl	L1 penalty for the factor loadings.
sigma0	The variance of the prior distribution of μ.
pi_alpha	The hyperparameter of the prior distribution of π. See details.
beta	A G by K0 matrix. Initial values for factor loadings (B). If null, beta will be initialized from normal distribution with mean zero and variance M0/K0. See details.
verbose	Whether to show some intermediate results. Default = False.
max_lambda	Whether to maximize over lambda. Default is True.
lambda	A M0 by K0 matrix. Initial values for the variances of factors. Each column is for a cell cluster. If null, lambda will initialize to be 1/M0. See details.
sigma	A G by M0 matrix. Initial values for the variance of idiosyncratic noises. Each column is for a cell cluster. If null, sigma will initialize to be 1. See details.
mu	A G by M0 matrix. Initial values for the gene expression mean of each cluster. Each column is for a cell cluster. If NULL, it will take the sample mean of cells weighted by the probability in each cluster. See details.
est_z	The iteration starts to update Z.
clus	Initial clustering of scRNASeq data. If NULL, the function will use PCA and Kmeans to do clustering initially.
p_min	Initialize parameters using genes expressed in at least p_min proportion of cells. If the number genes selected is less than min_gene, select min_gene genes with higest proportion of non zeros. Default = 0.4. If p_min is NULL, SIMPLE will estimate the dropout rate per gene and set 1-p_min to be the min(75% quantile of the dropout rate, 0.6)
cutoff	The value below cutoff is treated as no expression. Default = 0.1.
K	The number of PCs used in the initial clustering. Default is 20.
min_gene	Minimal number of genes used in the initial phase. Default: 2000. See details.
num_mc	The number of Gibbs steps for generating new imputed data after the parameters have been updated during Monte Carlo EM. Default = 3.
fix_num	If true, always use min_gene number of genes with the highest proportion of non zeros in the initial phase. Default = F. See details.
mcmc	The number of Gibbs steps to sample imputed data after EM. Default = 50.
burnin	The number of burnin steps before sample imputed data after EM. Default = 2.

Details

We assume that the cells come from M0 clusters. Within each cell cluster, the 'true' gene expression is modeled by a multivariate Gaussian distribution whose covariance matrix can be composed into a low rank matrix (a couple of latent gene modules) and idiosyncratic noises. Gene modules are shared among cell clusters, although the coexpression level of each gene module in different cell cluster can be different.
Suppose there are G genes and n cells. For each cell cluster, the gene expression follows Y|Z=m~MVN(μ_m, BΛ_m B^T + Σ_m) where B is a G by K0 matrix, Σ_m is a G by G diagonal matrix whose diagonal entries are specified by sigma, and Λ_m is a K0 by K0 diagonal matrix whose diagonal entries are specified by lambda. P(Z_m) = π_m where π~Dir(α).

The algorithm first runs Monte Carlo EM using only the genes with low dropout rate (initial phase) and initializes factor loadings and clustering membership. Then it runs more rounds of Monte Carlo EM using all the genes. In the initial phase, we use the genes with dropout rate less than 1 - p_min; if the number of such genes is less than min_gene, we rank the genes by the number cells with nonzero expression and keep the top min_gene genes. If fix_num is true, then we always keep the top min_gene genes in the initial phase.

After Monte Carlo EM, we obtain MLE of B, Λ_m, etc. If num_mc > 0, this function will sample multiple imputed values and cell factors at the MLEs of all the parameters by Gibbs sampling. Based on multiple imputed values, it will evaluate cluster stability for each cell (consensus_cluster). It will also output the posterior mean and variance for the imputed values and cell factors.

Value

SIMPLE returns a list of results in the following order.

1. loglik The log-likelihood of the full imputed gene expression at each iteration.

2. loglik_tot The log-likelihood of the full imputed gene expression at each iteration and the prior of B matrix.

3. BIC BIC which is -2 *loglik_tot + penalty on the number of parameters. Can be used to select paramters.

4. pi The prior probabilites of cells belong to each cluster.

5. mu Mean expression for each gene in each cluster

6. sigma Variances of idiosyncratic noises for each gene in each cluster.

7. beta Factor loadings.

8. lambda Variances of factors for each cluster.

9. z The posterior probability of each cell belonging to each cluster.

10. Yimp0 A matrix contains the expectation of gene expression specified by the model.

11. pg A G by M0 matrix, dropout rate for each gene in each cluster estimated from initial clustering.

12. initclus Output initial cluster results.

13. impt A matrix contains the mean of each imputed entry by sampling multiple imputed values while the parameters are MLE. If mcmc <= 0, output the imputed expressoin matrix at last step of EM

14. impt_var A matrix contains the variance of each imputed entry by sampling multiple imputed values while the parameters are MLE. NULL if mcmc <= 0.

15. Ef If mcmc >0, output posterior means of factors given observed data (a n by K0 matrix). If mcmc <= 0, output conditional expectation of the factors for each cluster E(f_i|z_i= m) at the last step of EM. A list with length M0, each element in the list is a n by K0 matrix.

16. Varf If mcmc >0, output posterior variances of factors given observed data (a n by K0 matrix). If mcmc <= 0, output conditional covariance matrix of factors for each cluster Var(f_i|z_i = m) at the last step of EM. A list with length M0, each element in the list is a K0 by K0 matrix.

17. consensus_cluster Score for the clustering stability of each cell by multiple imputations. NULL if mcmc <=0.

Author(s)

Zhirui Hu, [e-mail]

Songpeng Zu, [e-mail]

See Also

SIMPLE_B

Examples

library(foreach)
library(doParallel)
library(SIMPLEs)

# simulate number of clusters
M0 <- 3
# number of cells
n <- 300
# simulation_bulk and getCluster is defined in the util.R under the util directory of the corresponding github repository.
source("utils/utils.R")
simu_data <- simulation_bulk(n = 300, S0 = 20, K = 6, MC = M0, block_size = 32, indepG = 1000 - 32 * 6, verbose = F, overlap = 0)
Y2 <- simu_data$Y2
# number of factors
K0 <- 6
# parallel
registerDoParallel(cores = 6)
# estimate the parameters and sample imputed values
result <- SIMPLE(Y2, K0, M0, clus = NULL, K = 20, p_min = 0.5, max_lambda = T, min_gene = 200, cutoff = 0.01)
# evaluate cluster performance
celltype_true <- simu_data$Z
mclust::adjustedRandIndex(apply(result$z, 1, which.max), celltype_true)
# or redo clustering based on imputed values (sometimes work better for real data)
getCluster(result$impt, celltype_true, Ks = 20, M0 = M0)[[1]]

JunLiuLab/SIMPLEs documentation built on March 18, 2021, 3:10 a.m.