Contents

1 Introduction

Single-cell mRNA sequencing can uncover novel cell-to-cell heterogeneity in gene expression levels within seemingly homogeneous populations of cells. However, these experiments are prone to high levels of technical noise, creating new challenges for identifying genes that show genuine heterogeneous expression within the group of cells under study.

BASiCS (Bayesian Analysis of Single-Cell Sequencing data) is an integrated Bayesian hierarchical model that propagates statistical uncertainty by simultaneously performing data normalisation (global scaling), technical noise quantification and two types of supervised downstream analyses:

  1. For a single group of cells (Vallejos, Marioni, and Richardson 2015): BASiCS provides a criterion to identify highly (and lowly) variable genes within the group.
  2. For two (or more) groups of cells: BASiCS allows the identification of changes in gene expression between the groups. As in traditional differential expression tools, BASiCS can uncover changes in mean expression between the groups. Besides this, BASiCS can also uncover changes in gene expression variability in terms of:
  1. Over-dispersion (Vallejos, Richardson, and Marioni 2016) — a measure for the excess of cell-to-cell variability that is observed with respect to Poisson sampling, after accounting for technical noise. This feature has led, for example, to novel insights in the context of immune cells across aging (Martinez-Jimenez et al. 2017). However, due to the strong mean/over-dispersion confounding that is typically observed for scRNA-seq datasets, the assessment of changes in over-dispersion is restricted to genes for which mean expression does not change between the groups.
  2. Residual over-dispersion (Eling et al. 2018) — a residual measure of variability given by departures with respect to a global mean/over-dispersion trend. Positive residual over-dispersion indicates that a gene exhibits more variation than expected relative to genes with similar expression levels; negative residual over-dispersion suggests less variation than expected. To use this feature, please set Regression = TRUE as a function parameter in BASiCS_MCMC.

In all cases, a probabilistic output is provided and a decision rule is calibrated using the expected false discovery rate (EFDR) (Newton et al. 2004).

A brief description for the statistical model implemented in BASiCS is provided in Section 6 of this document. The original implementation of BASiCS (Vallejos, Marioni, and Richardson 2015) requires the use of spike-in molecules — that are artificially introduced to each cell’s lysate — to perform these analyses. More recently, Eling et al. (2018) extendeded BASiCS to also address datasets for which spikes-ins are not available (see Section 4). To use this feature, please set WithSpikes = FALSE as a function parameter in BASiCS_MCMC.

Important: BASiCS has been designed in the context of supervised experiments where the groups of cells (e.g. experimental conditions, cell types) under study are known a priori (e.g. case-control studies). Therefore, we DO NOT advise the use of BASiCS in unsupervised settings where the aim is to uncover sub-populations of cells through clustering.


2 Quick start

Parameter estimation is performed using the BASiCS_MCMC function. Downstream analyses implemented in BASiCS rely on appropriate post-processing of the output returned by BASiCS_MCMC. Essential parameters for running BASiCS_MCMC are:

If the optional parameter PrintProgress is set to TRUE, the R console will display the progress of the MCMC algorithm. For other optional parameters refer to help(BASiCS_MCMC).

Here, we illustrate the usage of BASiCS_MCMC using a built-in synthetic dataset.

Data <- makeExampleBASiCS_Data()
Chain <- BASiCS_MCMC(Data = Data, N = 1000, Thin = 10, Burn = 500, 
                     PrintProgress = FALSE, Regression = TRUE)
## -----------------------------------------------------
## MCMC sampler has been started: 1000 iterations to go.
## -----------------------------------------------------
## -----------------------------------------------------
## End of Burn-in period.
## -----------------------------------------------------
##  
## -----------------------------------------------------
## -----------------------------------------------------
## All 1000 MCMC iterations have been completed.
## -----------------------------------------------------
## -----------------------------------------------------
##  
## -----------------------------------------------------
## Please see below a summary of the overall acceptance rates.
## -----------------------------------------------------
##  
## Minimum acceptance rate among mu[i]'s: 0.424
## Average acceptance rate among mu[i]'s: 0.63544
## Maximum acceptance rate among mu[i]'s: 0.788
##  
##  
## Minimum acceptance rate among delta[i]'s: 0.438
## Average acceptance rate among delta[i]'s: 0.51916
## Maximum acceptance rate among delta[i]'s: 0.614
##  
##  
## Acceptance rate for phi (joint): 0.876
##  
##  
## Minimum acceptance rate among nu[j]'s: 0.424
## Average acceptance rate among nu[j]'s: 0.523133
## Maximum acceptance rate among nu[j]'s: 0.718
##  
##  
## Minimum acceptance rate among theta[k]'s: 0.704
## Average acceptance rate among theta[k]'s: 0.704
## Maximum acceptance rate among theta[k]'s: 0.704
##  
## -----------------------------------------------------
## 

As a default, the code above runs the original implementation mode of BASiCS (spikes without regression; see Section 4). To use the regression BASiCS model (Eling et al. 2018), please set Regression = TRUE. To use the no-spikes implementation of BASiCS, please add WithSpikes = FALSE as an additional parameter.

The Data and Chain (a BASiCS_Chain object) objects created by the code above can be use for subsequent downstream analyses. See Section 3.2 for highly/lowly variable gene detection (single group of cells, see also functions BASiCS_DetectHVG and BASiCS_DetectLVG) and Section 3.3 for differential expression analyses (two groups of cells, see also function BASiCS_TestDE).

Important remarks:

Typically, setting N=20000, Thin=20 and Burn=10000 leads to stable results.

3 Complete workflow

3.1 The input dataset

The input dataset for BASiCS must be stored as an SingleCellExperiment object (see SingleCellExperiment package).

The generation of the input SingleCellExperiment object requires a matrix of raw counts Counts (columns: cells, rows: genes) after quality control (e.g. removing low quality cells) and filtering of lowly expressed genes. If spike-in molecules are contained in Counts, a logical vector Tech is required to indicate which rows contain technical spike-in molecules and a data.frame object SpikeInfo containing the names of the spike-in molecules in the first column and the absolute number of molecules per well in the second column. More details are provided in section 3.1. If spike-ins are not available, a vector BatchInfo containing batch information is required.

3.1.1 With spike-in genes

The newBASiCS_Data function can be used to create the input data object based on the following information:

  • Counts: a matrix of raw expression counts with dimensions \(q\) times \(n\). Within this matrix, \(q_0\) rows must correspond to biological genes and \(q-q_0\) rows must correspond to technical spike-in genes. Gene names must be stored as rownames(Counts).

  • Tech: a logical vector (TRUE/FALSE) with \(q\) elements. If Tech[i] = FALSE the gene i is biological; otherwise the gene is spike-in. This vector must be specified in the same order of genes as in the Counts matrix.

  • SpikeInfo (optional): a data.frame with \(q-q_0\) rows. First column must contain the names associated to the spike-in genes (as in rownames(Counts)). Second column must contain the input number of molecules for the spike-in genes (amount per cell). If a value for this parameter is not provided when calling newBASiCS_Data, SpikeInfo is set as NULL as a default value. In those cases, the BatchInfo argument has to be provided to allow using the no-spikes implementation of BASiCS.

  • BatchInfo (optional): vector of length \(n\) to indicate batch structure (whenever cells have been processed using multiple batches). If a value for this parameter is not provided when calling newBASiCS_Data, BASiCS will assume the data contains a single batch of samples.

For example, the following code generates a synthetic dataset with 50 genes (40 biological and 10 spike-in) and 40 cells.

set.seed(1)
Counts <- matrix(rpois(50*40, 2), ncol = 40)
rownames(Counts) <- c(paste0("Gene", 1:40), paste0("Spike", 1:10))
Tech <- c(rep(FALSE,40),rep(TRUE,10))
set.seed(2)
SpikeInput <- rgamma(10,1,1)
SpikeInfo <- data.frame("SpikeID" = paste0("Spike", 1:10), 
                        "SpikeInput" = SpikeInput)

# No batch structure
DataExample <- newBASiCS_Data(Counts, Tech, SpikeInfo)

# With batch structure
DataExample <- newBASiCS_Data(Counts, Tech, SpikeInfo, 
                              BatchInfo = rep(c(1,2), each = 20)) 

To convert an existing SingleCellExperiment object (Data) into one that can be used within BASiCS, meta-information must be stored in the object.

  • isSpike(Data, "ERCC") <- Tech: the logical vector indicating biological/technical genes (see above) must be stored in the int_metadata slot via the isSpike function.

  • metadata(Data): the SpikeInfo object is stored in the metadata slot of the SingleCellExperiment object: metadata(Data) <- list(SpikeInput = SpikeInfo[,2], BatchInfo = BatchInfo).

  • colData(Data)$BatchInfo: the BatchInfo object is stored in the colData slot of the SingleCellExperiment object.

Once the additional information is included, the object can be used within BASiCS.

3.1.2 Without spike-in genes

To run BASiCS without incorporating reads from technical spike-in genes, and existing SingleCellExperiment object can be used. The only modification to the existing object is to assign the colData(Data)$BatchInfo slot.

set.seed(1)
CountsNoSpikes <- matrix(rpois(50*40, 2), ncol = 40)
rownames(CountsNoSpikes) <- paste0("Gene", 1:50)

# With batch structure
DataExampleNoSpikes <- SingleCellExperiment(assays = list(counts = CountsNoSpikes), 
                      colData = data.frame(BatchInfo = rep(c(1,2), each = 20))) 

Note: BASiCS assumes that a pre-processing quality control step has been applied to remove cells with poor quality data and/or lowly expressed genes that were undetected through sequencing. When analysing multiple groups of cells, the gene filtering step must be jointly applied across all groups to ensure the same genes are retained.

The function BASiCS_Filter can be used to perform this task. For examples, refer to help(BASiCS_Filter). Moreover, the scater package provides enhanced functionality for the pre-processing of scRNA-seq datasets.

3.2 Analysis for a single group of cells

We illustrate this analysis using a small extract from the MCMC chain obtained in (Vallejos, Richardson, and Marioni 2016) when analysing the single cell samples provided in (Grün, Kester, and Oudenaarden 2014). This is included within BASiCS as the ChainSC dataset.

data(ChainSC)

The following code is used to identify highly variable genes (HVG) and lowly variable genes (LVG) among these cells. The VarThreshold parameter sets a lower threshold for the proportion of variability that is assigned to the biological component (Sigma). In the examples below:

  • HVG are defined as those genes for which at least 60% of their total variability is attributed to the biological variability component.
  • LVG are defined as those genes for which at most 40% of their total variability is attributed to the biological variability component.

For each gene, these functions return posterior probabilities as a measure of HVG/LVG evidence. A cut-off value for these posterior probabilities is set by controlling the EFDR (as a default option, EFDR is set as 0.10).

par(mfrow = c(2,2))
HVG <- BASiCS_DetectHVG(ChainSC, VarThreshold = 0.6, Plot = TRUE)
LVG <- BASiCS_DetectLVG(ChainSC, VarThreshold = 0.2, Plot = TRUE)

To access the results of these tests, please use.

head(HVG$Table)
##     GeneIndex GeneName       Mu    Delta     Sigma      Prob   HVG
## 286       286   Rhox13 9.138169 2.417712 0.7706052 0.9733333  TRUE
## 250       250   Phlda2 9.372955 2.171966 0.7462231 0.9600000  TRUE
## 21         21    Amacr 7.164534 1.598386 0.6652893 0.6933333  TRUE
## 320       320    Smoc1 5.623408 1.805558 0.6700561 0.6800000 FALSE
## 122       122   Engase 7.082273 1.579305 0.6445464 0.6533333 FALSE
## 66         66   Cep170 3.712116 1.580848 0.6261821 0.5600000 FALSE
head(LVG$Table)
##     GeneIndex GeneName        Mu      Delta      Sigma Prob  LVG
## 20         20    Ahsa1  167.4669 0.06774259 0.09381064    1 TRUE
## 37         37   Atp5g2  344.9518 0.06381609 0.09037440    1 TRUE
## 55         55     Btf3  263.2987 0.04880375 0.06861493    1 TRUE
## 69         69   Chchd2  576.1569 0.03078477 0.04703723    1 TRUE
## 141       141    Fkbp4  334.4190 0.07084192 0.09782010    1 TRUE
## 147       147     Ftl1 2296.2110 0.04504962 0.06656222    1 TRUE

Note: this decision rule implemented in this function has changed with respect to the original release of BASiCS (where EviThreshold was defined such that EFDR = EFNR). However, the new choice is more stable (sometimes, it was not posible to find a threshold such that EFDR = EFNR).

3.3 Analysis for two groups of cells

To illustrate the use of the differential mean expression and differential over-dispersion tests between two cell populations, we use extracts from the MCMC chains obtained in (Vallejos, Richardson, and Marioni 2016) when analysing the (Grün, Kester, and Oudenaarden 2014) dataset (single cells vs pool-and-split samples). These were obtained by independently running the BASiCS_MCMC function for each group of cells.

data(ChainSC)
data(ChainRNA)
Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA,
                      GroupLabel1 = "SC", GroupLabel2 = "PaS",
                      EpsilonM = log2(1.5), EpsilonD = log2(1.5),
                      EFDR_M = 0.10, EFDR_D = 0.10,
                      Offset = TRUE, PlotOffset = FALSE, Plot = TRUE)

In BASiCS_TestDE, EpsilonM sets the log2 fold change (log2FC) in expression (\(\mu\)) and EpsilonD the log2FC in over-dispersion (\(\delta\)). As a default option: EpsilonM = EpsilonD = log2(1.5) (i.e. the test is set to detect absolute increases of 50% or above). To adjust for differences in overall mRNA content, an internal offset correction is performed when OffSet=TRUE. This is the recommended default setting.

The resulting output list can be displayed using

head(Test$TableMean)
##        GeneName MeanOverall   Mean1   Mean2 MeanFC MeanLog2FC ProbDiffMean
## 1 1700094D03Rik      58.539  54.352  62.726  0.854     -0.228        0.000
## 2 1700097N02Rik      37.421  36.919  37.922  0.946     -0.081        0.013
## 3 1810026B05Rik      14.738  10.459  19.018  0.558     -0.842        0.867
## 4 2310008H04Rik      18.844  15.267  22.422  0.672     -0.573        0.493
## 5 2410137M14Rik      19.050  16.674  21.426  0.768     -0.380        0.280
## 6 4930402H24Rik     962.255 967.739 956.771  1.004      0.006        0.000
##   ResultDiffMean
## 1         NoDiff
## 2         NoDiff
## 3           PaS+
## 4         NoDiff
## 5         NoDiff
## 6         NoDiff
head(Test$TableDisp)
##        GeneName MeanOverall DispOverall Disp1 Disp2 DispFC DispLog2FC
## 1 1700094D03Rik      58.539       0.225 0.321 0.128  2.307      1.206
## 2 1700097N02Rik      37.421       0.338 0.476 0.200  2.342      1.228
## 3 1810026B05Rik      14.738       0.359 0.434 0.284  1.378      0.463
## 4 2310008H04Rik      18.844       0.392 0.441 0.343  1.204      0.268
## 5 2410137M14Rik      19.050       0.396 0.491 0.302  1.663      0.734
## 6 4930402H24Rik     962.255       0.103 0.186 0.020  9.209      3.203
##   ProbDiffDisp      ResultDiffDisp
## 1        0.880                 SC+
## 2        0.813                 SC+
## 3        0.520 ExcludedFromTesting
## 4        0.453              NoDiff
## 5        0.640              NoDiff
## 6        1.000                 SC+

Due to the confounding between mean and over-dispersion that is typically observed in scRNA-seq datasets, the non-regression BASiCS model (run using Regression = FALSE as a function parameter in BASiCS_MCMC) can only be used to assess changes in over-dispersion for those genes in which the mean expression does not change between the groups. In this case, we recommend users to use EpsilonM = 0 as a conservative option to avoid changes in over-dispersion to be confounded by mean expression (the genes for which mean expression changes are marked as ExcludedFromTesting in the Test$TableDisp$ResultDiffDisp slot).

Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA,
                      GroupLabel1 = "SC", GroupLabel2 = "PaS",
                      EpsilonM = 0, EpsilonD = log2(1.5),
                      EFDR_M = 0.10, EFDR_D = 0.10,
                      Offset = TRUE, PlotOffset = FALSE, Plot = FALSE)

Note: If the regression BASiCS model has been used (Regression = TRUE as a function parameter in BASiCS_MCMC), BASiCS_TestDE will also report changes in residual over-dispersion (not confounded by mean expression) between the groups (see Section 4 in this vignette).

4 Alternative implementation modes

Beyond its original implementation, BASiCS has been extended to adress experimental designs in which spike-in molecules are not available as well as to address the confounding that is typically observed between mean and over-dispersion for scRNA-seq datasets (Eling et al. 2018). Alternative implementation modes are summarised below: