vignette.qmd

---
title: "Expression density estimation with Exprest"
author: "Simon Anders"
format: html
---

### The Task

In single-cell RNA-Seq data analysis, we are often interested in the distribution of a gene's expression strength over a population of cells or in comparing the expression distributions between several population. Usually, this is done with *violin plots*, as generated, for example, with Seurat's `VlnPlot` function. This should allow us to see whether expression is tightly controlled and hence approximately the same in all cells or whether it varies widely from cell to cell, and whether it is mono- or bimodal.

Usually, these violin plots are generated by performing a kernel density estimation (KDE) over log-normalized expression data, as performed by R's standard `density` function. However, the estimated distribution will be overly wide due to Poisson broadening, and the many zeroes in the data often give rise to a second mode at zero, even if the true underlying distribution is monomodal.

### Demonstration with a simple toy example

We use the generate a toy data set to illustrate how the standard approach leads to artifacts. 

In a typical microfluidic scRNA-Seq experiments, the total UMI count per cell (sum over all genes) is a few thousand. We simulate $n=1000$ cells and assign to each cell a total UMI count `s` drawn ar random such that the values of `s` are scattered around 5000.

```{r}
set.seed( 1234 )

n <- 1000
s <- round( 10^rnorm( n, mean=log10(5000), sd=.2 ) )

head( s, 10 )
```

We first consider a hypothetical gene G that has exactly the same expression in all the cells:
in every cell, 0.03% of all its mRNA molecules stem from this gene. For any given UMI count, the probablity that it is mapped to gene G is therefore $q=0.03\%=3\cdot10^4$. If a cell has `s` UMI counts, we can therefore draw from a binomial distribution with parameters `s` and `q` to get the counts `k`:

```{r}
k <- rbinom( n, s, 3e-4 )

head( k, 10 )
```

In a real RNA-Seq data set, `k` would be one row of the count matrix (corresponding to gene *G*) and `s` the column sums of the matrix.

We now perform log-normalization, using the formula $y=\log( k/s\cdot 10^4+1)$, which is the method used by default by Seurat's `NormalizeData` and also used in most scRNASeq analyses:

```{r}
y <- log( k/s * 1e4 + 1 )

head(y)
```

Here is a histogram of the log-normalized expression values, with the value corresponding to the log-transformed expression fraction of 0.03% marked with a pink line:

```{r}
hist( y, breaks=30, col="lightblue" )
abline( v = log( 3e-4 * 1e4 + 1 ), col="pink" )
```

Note how the histogram seems to suggest that the expression level is bimodal: some of the cells express the gene with a strength fluctuating around 0.03%, other seem to not express the gene at all. In reality, however, there is no such difference: in *all* our cells, 0.03% of the mRNA molecules were from gene G. Whether of these 1 or 2 molecules end up to be sequenced, or none, is a matter of random chance. 

A kernel density estimate performed with `density` gives a maybe even more misleading picture:

```{r}
naive_density_estimate <- density( y )
plot( naive_density_estimate, main="expression of gene G", col="blue" )
abline( v = log( 3e-4 * 1e4 + 1 ), col="pink" )
```
This (rotated and mirrored) is also what you would see in a typical violin plot:

```{r}
dens <- density(y)
plot( c( naive_density_estimate$y, -rev(naive_density_estimate$y) ), 
      c( naive_density_estimate$x, rev(naive_density_estimate$x) ), 
   type="l", xlim=c(-2,2), xlab="", ylab="", col="blue" )
```
The `exprest` package allows you to get the an estimate of the expression density that shows the actual monomodal expression, with its single narrow peak at 0.03%:

```{r}
library( exprest )

exprest_density_estimate <- est_dens( k, s )

plot( exprest_density_estimate, type="l", xlab="log10(lambda)", col="red" )

lines( naive_density_estimate$x/log(10)-4, naive_density_estimate$y, col="blue", lty="dashed" )
abline( v=log10(3e-4), col="pink" )
axis(3, log(exp(0:5)-1)/log(10)-4, 0:5 )
```

In this plot, we have marked the true expression with a pink line at the log10 of $3\cdot 10{-4}$. Just for completeness, we 
have also added Seurat's expression scale on top (which is different, because it uses a pseudocount and a scaling by $10^4$. See below for more details.)

Before we explain how `exprest` works, we show two more examples: a more involved simulated example, and an application to real data.

### Second example

Exprest's `make_test_data` function produces a data frame with example data for 1500 cells:

```{r}
df <- make_test_data()

head(df)
```

Here, the column `s` is the total counts per cell (the count-matrix column sums). It looks as follows:

```{r}
hist( log10( df$s ), breaks=30 )
```

For the true expression strength, we use this time a bimodal distribution:

```{r}
hist( log10(df$true_lambda), breaks=30, col="pink" )
```

The resulting counts for the simulated genes get, again, a prominent extra peak at zero:

```{r}
df$y <- log10( df$k/df$s + 1e-4 )
hist( df$y, breaks=30, freq=FALSE, col="lightblue" )
lines( density( df$y, bw=.1 ), col="blue" )
```

It is hard to see from this histogram (or from a corresponding violin plot of similar shape) that the gene's expression is actually bimodal and not trimodal.


The `exprest` estimate, however, resolves the issue:

```{r}
dens <- est_dens( df$k, df$s )

hist( log10( df$true_lambda ), border="gray", col="pink", breaks=30, freq=FALSE )
lines( density( df$y, bw=.1 ), col="blue", lty="dashed" )
lines( dens, col="red" )
```

Here, the red line is the density estimate. It is plotted over a histogram of the true expression values. 

### Application to real data

```{r}
seu <- schard::h5ad2seurat('~/Blood_TSP1_30_version2d_10X_smartseq_scvi_Nov122024.h5ad')
```