title: "Visualization of gene expression with Nebulosa (in Seurat)"

author: "Jose Alquicira-Hernandez"

package: Nebulosa

output:

BiocStyle::html_document:

self_contained: yes

toc: true

toc_float: true

toc_depth: 2

code_folding: show

vignette: >

Due to the sparsity observed in single-cell data (e.g. RNA-seq, ATAC-seq), the

visualization of cell features (e.g. gene, peak) is frequently affected and

unclear, especially when it is overlaid with clustering to annotate cell

types. `Nebulosa` is an R package to visualize data from single cells based on

kernel density estimation. It aims to recover the signal from dropped-out

features by incorporating the similarity between cells allowing a "convolution"

of the cell features.

For this vignette, let's use `Nebulosa` with the `Seurat` package.

First, we'll do a brief/standard data processing.

```{r import_libraries}

Let's download a dataset of 3k PBMCs (available from 10X Genomics). This same

dataset is commonly used in Seurat vignettes. The code below will download,

store, and uncompress the data in a temporary directory.

```{r download_and_untar_file}

Then, we can read the gene expression matrix using the `Read10X` from `Seurat`

```{r read_data}

Let's create a Seurat object with features being expressed in at least 3 cells

and cells expressing at least 200 genes.

```{r create_seurat_object}

Remove outlier cells based on the number of genes being expressed in each cell

(below 2500 genes) and expression of mitochondrial genes (below 5%).

```{r qc}

Let's use `SCTransform` to stabilize the variance of the data by regressing out

the effect of the sequencing depth from each cell.

```{r norm, message=FALSE, warning=FALSE}

Once the data is normalized and scaled, we can run a

_Principal Component Analysis_ (PCA) first to reduce the dimensions of our

data from 26286 features to 50 principal components. To visualize the principal

components, we can run a

_Uniform Manifold Approximation and Projection for Dimension Reduction_ (UMAP)

using the first 30 principal components to obtain a two-dimentional space.

```{r dim_red, message=FALSE, warning=FALSE}

To assess cell similarity, let's cluster the data by constructing a

_Shared Nearest Neighbor _(SNN) _Graph_ using the first 30 principal components

and applying the _Louvain algorithm_.

```{r clustering, message=FALSE, warning=FALSE}

The main function from `Nebulosa` is the `plot_density`. For usability, it

resembles the `FeaturePlot` function from `Seurat`.

Let's plot the kernel density estimate for `CD4` as follows

```{r plot_cd4}

For comparison, let's also plot a standard scatterplot using `Seurat`

```{r cd4_comparison}

By smoothing the data, `Nebulosa` allows a better visualization of the global

expression of CD4 in myeloid and CD4+ T cells. Notice that the "random"

expression of CD4 in other areas of the plot is removed as the expression of

this gene is not supported by many cells in those areas. Furthermore, CD4+

cells appear to show considerable dropout rate.

Let's plot the expression of CD4 with `Nebulosa` next to the clustering results

```{r cd4_and_clustering}

We can now easily identify that clusters `0` and `2` correspond to CD4+ T cells

if we plot CD3D too.

```{r plot_cd3d}

Characterize cell populations usually relies in more than a single marker.

Nebulosa allows the visualization of the joint density of from multiple

features in a single plot.

Users familiarized with PBMC datasets may know that CD8+ CCR7+ cells usually

cluster next to CD4+ CCR7+ and separate from the rest of CD8+ cells. Let's aim

to identify Naive CD8+ T cells. To do so, we can just add another gene to the

vector containing the features to visualize.

```{r fig.height=10}

`Nebulosa` can return a *joint density* plot by multiplying the densities

from all query genes by using the `joint = TRUE` parameter:

```{r fig.height=14}

When compared to the clustering results, we can easily identify that Naive

CD8+ T cells correspond to cluster `8`.

`Nebulosa` returns the density estimates for each gene along with the joint

density across all provided genes. By setting `combine = FALSE`, we can obtain

a list of ggplot objects where the last plot corresponds to the joint density

estimate.

Likewise, the identification of Naive CD4+ T cells becomes straightforward by

combining `CD4` and `CCR7`:

```{r fig.height=14}

Notice that these cells are mainly constrained to cluster `0`

```{r fig.height=10}

In summary,`Nebulosa`can be useful to recover the signal from dropped-out genes

and improve their visualization in a two-dimensional space. We recommend using

`Nebulosa` particularly for dropped-out genes. For fairly well-expressed genes,

the direct visualization of the gene expression may be preferable. We encourage

users to use `Nebulosa` along with the core visualization methods from the

`Seurat` and `Bioconductor` environments as well as other visualization methods

to draw more informed conclusions about their data.