Estimation of capture efficiencies

Apart from the raw data, another parameter which user may need to provide is the mean capture efficiency  < β> and hence β can be further calculated using scaling factors estimated from any other methods. The default β is calculated to be total counts normalized to 6%. Or you can use BetaFun in bayNorm to estimate β:

data('EXAMPLE_DATA_list')
#Suppose the input data is a SummarizedExperiment object:
#Here we just used 30 cells for illustration
rse <- SummarizedExperiment::SummarizedExperiment(assays=SimpleList(counts=EXAMPLE_DATA_list$inputdata[,seq(1,30)]))
#SingleCellExperiment object can also be input in bayNorm:
#rse <- SingleCellExperiment::SingleCellExperiment(assays=list(counts=EXAMPLE_DATA_list$inputdata))

BETA_est<-BetaFun(Data=rse,MeanBETA=0.06)
summary(BETA_est$BETA)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.02299 0.03859 0.05255 0.06000 0.07473 0.11168

The function BetaFun selects a subset of genes such that outlier genes and genes with high proportion of zeros across cells are excluded. Then the total counts of the subset of genes in each cell are normalized to  < β>. Another way is to normalize scaling factors estimated from the R package scran to  < β>.

Run bayNorm

BETA_vec can be set to be NULL, then the β are estimated to be total counts normalized to 6%.

#Return 3D array normalzied data, draw 20 samples from posterior distribution:
bayNorm_3D<-bayNorm(
    Data=rse,
    BETA_vec = NULL,
    mode_version=FALSE,
    mean_version = FALSE,S=20
    ,verbose =FALSE,
    parallel = TRUE)

##   |                                                                              |                                                                      |   0%  |                                                                              |====                                                                  |   5%  |                                                                              |=======                                                               |  10%  |                                                                              |==========                                                            |  15%  |                                                                              |==============                                                        |  20%  |                                                                              |==================                                                    |  25%  |                                                                              |=====================                                                 |  30%  |                                                                              |========================                                              |  35%  |                                                                              |============================                                          |  40%  |                                                                              |================================                                      |  45%  |                                                                              |===================================                                   |  50%  |                                                                              |======================================                                |  55%  |                                                                              |==========================================                            |  60%  |                                                                              |==============================================                        |  65%  |                                                                              |=================================================                     |  70%  |                                                                              |====================================================                  |  75%  |                                                                              |========================================================              |  80%  |                                                                              |============================================================          |  85%  |                                                                              |===============================================================       |  90%  |                                                                              |==================================================================    |  95%  |                                                                              |======================================================================| 100%

#Return 2D matrix normalized data (MAP of posterior):
#Simply set mode_version=TRUE, but keep mean_version=FALSE


#Return 2D matrix normalized data (mean of posterior):
#Simply set mean_version=TRUE, but keep mode_version=FALSE

Non-UMI scRNAseq dataset

bayNorm’s mathematical model is suitable for UMI dataset. However it can be also applied on non-UMI dataset. In bayNorm, you need to specify the following parameter: * UMI_sffl: bayNorm can also be applied on the non-UMI dataset. However, user need to provide a scaled number. Raw data will be divided by the scaled number and bayNorm will be applied on the rounded scaled data. By doing so, the Dropout vs Mean expression plots will be similar to that of UMI dataset.

Generate 3D array or 2D matrix with existing estimated prior parameters.

If you have run bayNorm on a dataset before but want to output another kind of data (3D array or 2D matrix), you can use the function bayNorm_sup. It is important to input the existing estimated parameters by specifying the following parameter in bayNorm_sup:

• BETA_vec: If Conditions has been specified previously, then input unlist(bayNorm_output$BETA)
• PRIORS: input bayNorm_output$PRIORS_LIST
• Conditions: make sure to specify the same Conditions as before. You can find these two objects from the previous output of bayNorm function, which is a list.

#Now if you want to generate 2D matrix (MAP) using the same prior
#estimates as generated before:
bayNorm_2D<-bayNorm_sup(
    Data=rse,
    PRIORS=bayNorm_3D$PRIORS,
    input_params=bayNorm_3D$input_params,
    mode_version=TRUE,
    mean_version = FALSE,
    verbose =FALSE)

#Or you may want to generate 2D matrix 
#(mean of posterior) using the same prior
#estimates as generated before:
bayNorm_2D<-bayNorm_sup(
    Data=rse,
    PRIORS=bayNorm_3D$PRIORS,
    input_params=bayNorm_3D$input_params,
    mode_version=FALSE,
    mean_version = TRUE,
    verbose =FALSE)

Rationale of bayNorm

A scRNA-seq dataset is typically represented in a matrix of dimension P × Q, where P denotes the total number of genes observed and Q denotes the total number of cells studied. The element x_ij (i ∈ {1, 2, …, P} and j ∈ {1, 2, …, Q}) in the matrix represents the number of transcripts reported for the i^th gene in the j^th cell. This is equal to the total number of sequencing reads mapping to that gene in that cell for a non-UMI protocol. For UMI based protocols this is equal to the number of individual UMIs mapping to each gene. The matrix can include data from different groups or batches of cells, representing different biological conditions. This can be represented as a vector of labels for the cell groups or conditions (C_j). bayNorm generates for each gene (i) in each cell (j) a posterior distribution of original expression counts (x_ij⁰), given the observed scRNA-seq read count for that gene (x_ij).

A common approach for normalizing scRNA-seq data is based on the use of a global scaling factor (s_j), ignoring any gene specific biases. The normalized data x̃_ij is obtained by dividing the raw data for each cell j by the its global scaling factor s_j:

In bayNorm, we implement global scaling using a Bayesian approach. We assume given the original number of transcripts in the cell (x_ij⁰), the number of transcripts observed (x_ij) follows a Binomial model with probability β_j, which we refer to as capture effeiciency and it represents the probability of original transcripts in the cell to be observed. In addition, we assume that the original number or true count of the i^th gene in the j^th cell (x_ij⁰) follows Negative Binomial distribution with parameters mean (μ), size (or dispersion parameter, ϕ), such that: $$\Pr(x^0_{ij}=n|\phi_i,\mu_i)=\frac{\Gamma(n+\phi_i)}{\Gamma(\phi_i)n!}(\frac{\phi_i}{\mu_i+\phi_i})^{\phi_i}(\frac{\mu_i}{\mu_i+\phi_i})^{n}$$

So, overall we have the following model:

Using the Bayes rule, we have the following posterior distribution of original number of mRNAs for each gene in each cell:

The prior parameters μ and ϕ of each gene were estimated using an empirical Bayesian method as discussed in detail in Supplementary Information of 1.

For more details about the rationale of bayNorm, please see 1.

Estimation of capture efficiencies

Cell specific capture efficiency β_j and global scaling factor (s_j) are closely related. We can transform scaling factors estimated by different methods (see below) into β_j values with the following formula:

β̄ or  < β>, a scalar, is an estimate of global mean capture efficiency across all cells, which ranges between 0 and 1.

There are two different methods for estimating β̄ and β_j:

1. If spike-ins or smFISH data are available they can be used to estimate capture efficiencies. We can either divide the total number of observed spik-ins in each cell by the total number of input spike-ins, or we can fit a linear regression to estimate the cell specific β_j. If smFISH data is available, we can fit a linear regression between the mean expression of raw data (response variable) and the mean expression of the smFISH data (explanatory variable). The coefficient of the explanatory variable can be used as β̄.
2. The raw data itself can be directly used for estimation of cell specific global scaling factors (s_j). Then an estimate of β̄ can be used to estimate β_j. There are different methods available for estimation of global scaling factors. Some were developed for bulk RNA-seq data and some are specific to scRNA-seq data. The value of β̄ depends on the protocol used and can be batch dependent. For example, for Droplet based protocol, it is about 0.06 (inDrop) or 0.12 (Drop-seq). β̄ can also be estimated by spike-ins or smFISH data as explained above.

We finally note, that estimates of capture efficiency discussed above will assume cells have simular original transcript content. Therefore, the bayNorm outputs estimates of original transcript counts for a typical cell, which is corrected for variation in cell size and transcript content. This is usually desirable for down-stream analysis such as DE detection. However, if one is interested in absolute origianl count and has additional information such as cell size or total transcirpt content per cell, the capture efficiencies can be approporiatly rescaled for this purpose.

Estimation of prior parameters

Method of Moments

A heuristic way of estimating μ_i and ϕ_i is through a variant of the Method of Moments. The first step is to do a simple normalization of the raw data, to scale expressions given the cell specific capture efficiencies (β_j). The simple normalized count x_ij^s is calculated as following:

where the numerator of the scaling factor of x_ij is obtained by taking the average of scaled total counts across cells.

Based on simple normalized data, we are able to estimate prior parameters μ and ϕ of the Negative Binomial distribution using the Method of Moments Estimation (MME), which simply equates the theoretical and empirical moments. This estimation method is fast and simulations suggests it provides good estimates of μ but the drawback is that the estimation of ϕ show a systematic bias (see Supplementary Figure S24 a-b in 1).

The combined method (default setting in bayNorm for estimating priors)

Based on simulation studies (Supplementary Figure S24 in 1), the most robust and efficient estimation of μ and ϕ can be obtained using the following combined approach, which is the default setting in bayNorm:

1. Based on simple normalized data, we use the MME method for each gene to obtain MME estimated μ and ϕ.
2. Although the BB estimated ϕ is much closer to the true ϕ, many estimates are at the upper boundary of the search space (Supplementary Figures S24 c-d in 1). So, we find adjusting the MME estimated ϕ by a factor which can be estimated by fitting a linear regression between MME estimated ϕ and BB estimated ϕ works best (Supplementary Figures S24 c-d in 1). This adjusted MME estimated ϕ together with the MME estimated μ and estimates of β_j can be used in approximating posterior distribution for each gene in each cell.

Cells are grouped together for prior estimation, based on cell-specific attributes (C_j). Prior estimation can be done over all cells irrespective of the experimental condition. We refer to this procedure as “global”. Alternatively, suppose that there are multiple groups of cells in the datasets and we have reasons to believe each group could behave differently. Then we can estimate the prior parameters “μ and ϕ” within each group respectively (within groups with the same C_j value). We refer to this procedure as “local”. Estimating prior parameters across a certain group of cells based on “global” procedure allow for removing potential batch effects. Multiple groups normalization based on “local” procedure allows for amplifying the inter-groups’ differences while mitigating the intra-group’s variability, which is suitable for DE detection.