A UNIFIED STATISTICAL FRAMEWORK FOR SINGLE CELL AND BULK RNA SEQUENCING DATA

Abstract

Recent advances in technology have enabled the measurement of RNA levels for individual cells. Compared to traditional tissue-level bulk RNA-seq data, single cell sequencing yields valuable insights about gene expression profiles for different cell types, which is potentially critical for understanding many complex human diseases. However, developing quantitative tools for such data remains challenging because of high levels of technical noise, especially the "dropout" events. A "dropout" happens when the RNA for a gene fails to be amplified prior to sequencing, producing a "false" zero in the observed data. In this paper, we propose a Unified RNA-Sequencing Model (URSM) for both single cell and bulk RNA-seq data, formulated as a hierarchical model. URSM borrows the strength from both data sources and carefully models the dropouts in single cell data, leading to a more accurate estimation of cell type specific gene expression profile. In addition, URSM naturally provides inference on the dropout entries in single cell data that need to be imputed for downstream analyses, as well as the mixing proportions of different cell types in bulk samples. We adopt an empirical Bayes' approach, where parameters are estimated using the EM algorithm and approximate inference is obtained by Gibbs sampling. Simulation results illustrate that URSM outperforms existing approaches both in correcting for dropouts in single cell data, as well as in deconvolving bulk samples. We also demonstrate an application to gene expression data on fetal brains, where our model successfully imputes the dropout genes and reveals cell type specific expression patterns.

Keywords: EM algorithm; Gibbs sampling; Single cell RNA sequencing; empirical Bayes; hierarchical model.

Fig. 1

Plate model of URSM, with both single cell data (on the left) and bulk samples (on the right). The two greyed nodes X and Y represent observed gene expression levels. Node S is a binary variable representing dropout status in single cells, and node W represents the mixing proportions in bulk samples. The node π representing observation probability is double-circled because it is deterministic, and all model parameters are shown without circles, including the profile matrix A that links the two data sources.

Fig. 2

(a) Simulated logistic dropout probability curves for 100 single cells, as defined in equation (4.1). (b)–(d) True profile matrix A versus the estimated Â, plotted in the log scale, using (i) the naive sample mean estimation [equation (4.2)], (ii) a submodel using only single cell data and (iii) URSM with both single cell and bulk data. The L₁ loss Σ_{i, k} |Â_ik – A_ik| is reported on the top.

Fig. 3

The L₁ loss of recovering (a) the profile matrix, Σ_{i, k} |Â_ik – A_ik|, and (b) mixing proportions, Σ_{k, j} |Ŵ_kj – W_kj |. We evaluate DSA and ssNMF when the marker genes are extracted from single cell data using different thresholds of p-values, as well as under the oracle condition where the true marker genes are given. We evaluate Cibersort on estimating W when the input signature matrix is based on the estimated  from URSM. We report its performance when the entire  is used (“Cibersort all”), as well as when only the estimated marker genes are used (“Cibersort”). The performance of URSM is plotted with a square in both panels, which does not depend on thresholding p-values.

Fig. 4

(a) ROC curves of identifying dropout entries in single cell data. (b) True profile matrix A versus the sample average of a rank-3 NMF approximation, plotted in the log scale. The L₁ loss Σ_{i, k} |Â_ik – A_ik| is reported on the top.

Fig. 5

The average per cell type L₁ loss of recovering the profile matrix A and the mixing proportions W in 10 repetitions, with the standard deviations shown by the error bars, when (a) K_sc, K_bk ∈ {3, 4, 5} with N = 200 genes; (b) N = {200, 500, 1000} with K_sc = K_bk = 3. Each figure shows the performance of (i) URSM, (ii) DSA and ssNMF with 5 true marker genes and 3 imperfect marker genes per cell type as input and (iii) DSA and ssNMF under the oracle scenario where 5 true marker genes per cell type are provided. We also report the performance of Cibersort for estimating W using the estimated Â_unif from URSM as the input signature matrix.

Fig. 6

(a) Single cell gene expressions (log₂(FPKM + 1)) after removing 7 ambiguously labeled cells. Rows are 213 cells and columns are 273 genes. (b) PCA applied on the original single cell data with 220 labeled cells using 273 PC genes, where the Monocle algorithm is applied to construct pseudo developmental times. 7 cells are identified to be ambiguously labeled and are removed from our analyses (marked as triangles). (c) Entries in cleaned single cell data that are inferred to be dropout and imputed (marked in blue) versus the entries that are inferred to be structural zeros (marked in white) in cleaned single cell data. The entries with positive expression levels have no need for posterior inference, and are marked in grey. (d) After imputing dropout genes, PCA is conducted on the 213 cells using 273 PC genes, and the three different types of cells are more clearly separated.

Fig. 7

Deconvolution of bulk samples into three cell types, using (a) URSM, (b) Cibersort, (c) Digital Sorting Algorithm (DSA) and (d) semi-supervised Nonnegative Matrix Factorization (ssNMF).