---
title: "Transcriptional variability analysis of CD4^+^ T cells using BASiCS"
author: 
  - Alan O'Callaghan
  - Nils Eling
  - John C. Marioni
  - Catalina A. Vallejos
bibliography: Workflow.bib
urlcolor: Orange
output:
  bookdown::pdf_document2:
    fig_width: 6
    fig_height: 4
---


```{r setup_knitr, include = FALSE, cache = FALSE}
library("BiocStyle")
library("knitr")
knitr::opts_chunk$set(
  message = FALSE, error = FALSE, warning = FALSE,
  cache = 0, cache.path = "cache_supp/",
  fig.path = "figure_supp/"
)
```

# 1. Introduction

To date, three versions of the `r Biocpkg("BASiCS")` model have 
been proposed [@Vallejos2015BASiCS] [@Vallejos2016] [@Eling2018]. They differ in 
how inference is performed (e.g. using different priors) and the type of 
downstream analysis that is enabled by the model. 

- **Vallejos et al (2015)** [@Vallejos2015]: the original model uses 
information from extrinsic spike-in molecules (e.g. those introduced by 
[@Rna2005]) as *control features* to quantify technical noise. This enables the 
estimation of two sets of cell-specific normalisation parameters ($s_j$ and 
$\phi_j$) capturing technical (e.g. amplification biases) and biological 
(e.g. mRNA content) systematic differences across cells [@Vallejos2017]. A
probabilistic decision rule (based on $\delta_i$) was proposed to identify 
*highly variable genes* (HVGs) that capture the major sources of heterogeneity 
within the analysed cells [@Brennecke2013]. HVG detection is often used to 
perform feature selection, choosing the input set of genes for subsequent 
analyses. A similar rule was developed to highlight *lowly variable genes* 
(LVGs) that exhibit stable expression across the population of cells. These may 
relate to essential cellular functions and can assist the development of new 
data normalisation or integration strategies [@Lin2019].

- **Vallejos et al (2016)** [@Vallejos2016]: the model was extended to enable
**differential expression** analyses between two pre-specified groups of cells 
(e.g. different experimental conditions or cell types). This is achieved by 
comparing the posterior distribution of gene-specific parameters ($\mu_i$ and 
$\delta_i$). While several differential expression tools were previously 
proposed for scRNA-seq data (e.g. [@Kharchenko2014; @Finak2015]), some evidence 
suggests that these do not generally outperform popular bulk RNA-seq tools 
[@Soneson2018]. Moreover, most of these methods are only designed to uncover 
changes in overall expression, ignoring the more complex patterns that can arise 
at the single cell level [@Lahnemann2020]. Instead, `r Biocpkg("BASiCS")` 
embraces the high granularity of scRNA-seq data, uncovering changes in 
cell-to-cell transcriptional variability. As noted by @Vallejos2016, the inverse
relationship that is observed between mean expression and over-dispersion 
derived from (bulk and) scRNAseq can affect the interpretation of such analyses.
In particular, genes that are differentially expressed between two groups of
cells are likely to exhibit changes in both mean expression and variability,
due to the inverse relationship between these two quantities. Thus, 
comparisons of variability between populations must be restricted to genes
that do not exhibit changes in mean expression.

- **Eling et al (2018)** [@Eling2018]: the model was extended to account for
the strong relationship that is typically observed 
between gene-specific mean expression and over-dispersion estimates.
Eling *et al.*  [@Eling2018] introduced a *joint prior* specification for 
these parameters. This joint prior assumes that genes with similar mean 
expression ($\mu_i$) have similar over-dispersion parameters $\delta_i$. 
Effectively, this shrinks over-dispersion estimates towards a global trend
that captures the relationship between mean and over-dispersion (Figure XX). 
This improves posterior inference for over-dispersion parameters when the data 
is less informative (e.g. small sample size, lowly expressed genes) [@Eling2018].
This information-sharing approach is conceptually similar to that performed by 
@Love2014 and others, where sparse data is pooled to obtain more reliable 
estimates.
The global trend is then used to derive gene-specific *residual over-dispersion* 
parameters $\epsilon_i$ that are not confounded by mean expression. 
Similar to the DM values implemented in `r Biocpkg("scran")`, these are defined 
as deviations with respect to the overall trend (Figure XX).
`r Biocpkg("BASiCS")` also provides a probabilistic decision rule to 
perform differential expression analyses between two pre-specified 
groups of cells [@Vallejos2016; @Eling2018].
Furthermore, the model was extended using a horizontal integration framework
to allow its use in the absence of spike-in genes.
This is useful
for droplet-based scRNAseq protocols, given that it is not possible to ensure
that each droplet contains a specified quantity of spike-in molecules.
In this horizontal integration framework, technical variation is quantified
using replication [@Carroll2005]. In the absence of true technical replicates,
we assume that population-level characteristics of the cells are replicated
using appropriate experimental design. This requires that cells from the same
population have been randomly allocated to different batches. Given appropriate
experimental design, `r Biocpkg("BASiCS")` assumes that biological effects
are shared across batches, while technical variation leads to spurious
differences between cells in different batches.
It is this version of
the model that we focus on here, and that we recommend for most
users. Previous versions of the model are available within the package, but
are primarily useful for reproducibility purposes or for analysing
datasets that contain spike-in genes.


### Testing for changes in mean expression and over-dispersion

Differential mean and differential over-dispersion testing is done by computing
the tail posterior propabilities of the difference in mean expression or 
over-dispersion between two conditions ($A$ and $B$) being larger than an
evidence threshold $\tau_0$ or $\omega_0$ [@Bochkina2007, @Vallejos2016]:
$$
\begin{aligned}
  &\text{P}\left(\log\left(\frac{\mu_i^{(A)}}{\mu_i^{(B)}}\right) > \tau_0 | \text{Data}\right) > \alpha_m\\
  &\text{P}\left(\log\left(\frac{\delta_i^{(A)}}{\delta_i^{(B)}}\right) > \omega_0 | \text{Data}\right) > \alpha_d
\end{aligned}
$$

If the tail posterior probability is larger than a given propability threshold 
$\alpha_m$ or $\alpha_d$, the gene is considered to be differentially expressed
or differentially over-dispersed [@Vallejos2016]. The evidence threshold is
usually fixed *a priori* and the probability threshold is defined to
control the expected false discovery rate (EFDR) to 
(e.g. 10%) [@Newton2004, Vallejos2016].

As described by @Vallejos2016, estimates of the
over-dispersion parameters $\delta_i$ are negatively correlated to mean 
expression $\mu_i$. This indicates that in homogeneous populations of cells,
highly expressed genes tend to be less noisy than lowly expressed genes. 
Differential over-dispersion testing is therefore confounded by mean expression
changes. When assessing changes in over-dispersion $\delta_i$, only genes with
no changes in mean expression are considered (see @Vallejos2016).

### Correcting the mean-variability confounding effect

@Eling2018 extended BASiCS to account for the
confounding effect between mean expression and expression variability. 
For this purpose, we capture the relationship between mean and over-dispersion
parameters by introducing the following joint prior distribution for 
$(\mu_i, \delta_i)'$:
$$
  \mu_i \sim \text{log-Normal}\left(0, s^2_{\mu}\right), \hspace{0.5cm}
  \delta_i | \mu_i \sim \text{log-t}_{\eta}\left( \text{f}(\mu_i), \sigma^2 \right).
$$

The latter is equivalent to the non-linear regression model:
$$
  \log(\delta_i) = \text{f}(\mu_i) + \epsilon_i, \hspace{0.5cm} \epsilon_i \sim{}\text{t}_{\eta}(0,\sigma^2), 
$$

where $\mbox{f}(\mu_i)$ represents the over-dispersion (on the log-scale) that
is predicted by the global trend (across all genes) for a given mean expression 
$\mu_i$. Therefore, $\epsilon_i$ can be interpreted as a gene-specific 
*residual over-dispersion* parameter. Positive values for
$\epsilon_i$ indicate more variation than expected for genes with similar 
expression level. Similarly, negative values of $\epsilon_i$ suggest less 
variation than expected [@Eling2018].

In line with the probabilistic approach described above, we identified
statistically significant differences in residual over-dispersion for those 
genes where the tail posterior probability of observing a large difference 
between $\epsilon_i^A$ and $\epsilon_i^B$ exceeds a certain threshold:
$$
  \text{P}(\mid\epsilon_i^{A} - \epsilon_i^{B}\mid > \psi_0 \mid \text{Data}) > \alpha_R,
$$
where $\psi_0 > 0$ defines the minimum tolerance threshold. 
As a default choice, we assume $\psi_0 = \log_2(1.5) / \log_2(e) \approx 0.41$,
which translates into a 50\% increase in over-dispersion. 
The posterior probability threshold $\alpha_R$ is chosen to control the EFDR
(e.g. 10\%) [@Newton2004]. To support interpretability of the results,
we exclude genes that are not expressed in at least 2 cells per condition from
differential variability testing.
