---
title: "PsiNorm: a scalable normalization for single-cellRNA-seq data"
author: 
    name: Matteo Borella, Chiara Romualdi, and Davide Risso
    affiliation: Department of Biology and Department of Statistical Sciences, University of Padova 
date: "Last modified: April 13, 2021; Compiled: `r format(Sys.time(), '%B %d, %Y')`"
vignette: >
  %\VignetteEncoding{UTF-8}
output:
      BiocStyle::html_document  
---

<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{PsiNorm normalization}
-->

```{r, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


# Introduction

```{r, message=FALSE, warning=FALSE}
library(SingleCellExperiment)
library(splatter)
library(scater)
library(cluster)
library(scone)
```

PsiNorm is a scalable between-sample normalization for single cell RNA-seq count data based on the power-law Pareto type I distribution. It can be demonstrated that the Pareto parameter is inversely proportional to the sequencing depth, it is sample specific and its estimate can be obtained for each cell independently. PsiNorm computes the shape parameter for each cellular sample and then uses it as multiplicative size factor to normalize the data. The final goal of the transformation is to align the gene expression distribution especially for those genes characterised by high expression. Note that, similar to other global scaling methods, our method does not remove batch effects, which can be dealt with downstream tools.

To evaluate the ability of PsiNorm to remove technical bias and  reveal the true cell similarity structure, we used both an unsupervised and a supervised approach.
We first simulate a scRNA-seq experiment with four known clusters using the _splatter_ Bioconductor package.  Then in the unsupervised approach, we i) reduce dimentionality using PCA, ii) identify clusters using the _clara_ partitional method and then we iii)  computed the Adjusted Rand Index (ARI) to compare the known and the estimated partition.

In the supervised approach, we i) reduce dimentionality using PCA, and we  ii)  compute the silhouette index of the  known partition in the reduced dimensional space.

# Citation

If you use `PsiNorm` in publications, please cite the following article:

Borella, M., Martello, G., Risso, D., & Romualdi, C. (2021). PsiNorm: a scalable normalization for single-cell RNA-seq data. bioRxiv. https://doi.org/10.1101/2021.04.07.438822.

# Data Simulation

We simulate a matrix of counts with 2000 cellular samples and 10000 genes with splatter.

```{r}
set.seed(1234)
params <- newSplatParams()
N=2000
sce <- splatSimulateGroups(params, batchCells=N, lib.loc=12,
                           group.prob = rep(0.25,4),
                           de.prob = 0.2, de.facLoc = 0.06,
                           verbose = FALSE) 
```

`sce` is a SingleCellExperiment object with a single batch and four different cellular groups.

To visualize the data we used the first two Principal Components estimated starting from the raw log-count matrix.

```{r}
set.seed(1234)
assay(sce, "lograwcounts") <- log1p(counts(sce))
sce <- runPCA(sce, exprs_values="lograwcounts", scale=TRUE, ncomponents = 2)
plotPCA(sce, colour_by="Group")
```

# PsiNorm data normalization 

# Data Normalization with PsiNorm 

To normalize the raw counts we used the PsiNorm normalization and we visualized the data using the first two principal components. 

```{r}
sce<-PsiNorm(sce)
sce<-logNormCounts(sce)
head(sizeFactors(sce))
```

Note that running the `PsiNorm` function computes a set of size factors that are added to the SingleCellExperiment object.

The `logNormCounts` function can be then used to normalize the data by multiplying the raw counts and the size factors.

```{r}
set.seed(1234)
sce<-runPCA(sce, exprs_values="logcounts", scale=TRUE, name = "PsiNorm_PCA",
            ncomponents = 2)
plotReducedDim(sce, dimred = "PsiNorm_PCA", colour_by = "Group")
```

We can appreciate from the plot that PsiNorm allows a better separation among known cellular groups.

## Unsupervised approach: Adusted Rand Index  

We calculate ARI of both raw counts and PsiNorm normalized counts after PCA dimension reduction and $clara$ clustering (with $k$ equal to the simulated number of clusters); higher the ARI, better the normalization. 

```{r}
groups<-cluster::clara(reducedDim(sce, "PCA"), k=nlevels(sce$Group))
a<-paste("ARI from raw counts:", 
         round(mclust::adjustedRandIndex(groups$clustering, sce$Group), 
               digits = 3))

groups<-cluster::clara(reducedDim(sce, "PsiNorm_PCA"), k=nlevels(sce$Group))
b<-paste("ARI from PsiNorm normalized data:",
         round(mclust::adjustedRandIndex(groups$clustering, sce$Group), 
               digits = 3))

kableExtra::kable(rbind(a,b), row.names = FALSE)
```

Pareto normalization considerably increases the ARI index.

# Supervised approach: Silhouette index

We calculate the Silhouette index of both raw counts and PsiNorm normalized counts after tSNE dimension reduction exploiting known simulated clusters; higher the Silhouette, better the normalization. 

```{r}
dist<-daisy(reducedDim(sce, "PCA"))
dist<-as.matrix(dist)
a<-paste("Silhouette from raw counts:", round(summary(
    silhouette(x=as.numeric(as.factor(sce$Group)),
               dmatrix = dist))$avg.width, digits = 3))

dist<-daisy(reducedDim(sce, "PsiNorm_PCA"))
dist<-as.matrix(dist)
b<-paste("Silhouette from PsiNorm normalized data:", round(summary(
    silhouette(x=as.numeric(as.factor(sce$Group)),
               dmatrix = dist))$avg.width, digits = 3))
kableExtra::kable(rbind(a,b), row.names = FALSE)
```

Pareto normalization considerably increases the Silhouette index. 

# Correlation of PC1 and PC2 with sequencing depth

To check if PsiNorm is able to capture technical noise and remove unwanted variation within a dataset (due for instance to differences in sequencing depth), we check whether the first two PCs are capturing technical variance. We computed the maximum correlation obtained between PC1 and PC2 and cell sequencing depths; a higher correlation indicates that the normalization was not able to properly remove noise.

```{r}
set.seed(4444)
PCA<-reducedDim(sce, "PCA") 
PCAp<-reducedDim(sce, "PsiNorm_PCA")
depth<-apply(counts(sce), 2, sum)
a<-paste("The Correlation with the raw data is:",
            round(abs(max(cor(PCA[,1], depth), cor(PCA[,2], depth))), digits=3))
b<-paste("The Correlation with the PsiNorm normalized data is:",
            round(abs(max(cor(PCAp[,1], depth), cor(PCAp[,2], depth))), digits = 3))
kableExtra::kable(rbind(a,b), row.names = FALSE)
```

Our results demonstrate that the correlation significantly decreases after the PsiNorm normalization.

# Using PsiNorm in `scone()`

As for other normalizations, `scone` includes a wrapper function to use PsiNorm in the SCONE evaluation framework.
See Section 3.2 of the "Introduction to SCONE" vignette for an example on how to use PsiNorm within the main `scone()` function.

# Using PsiNorm with Seurat

The PsiNorm normalization method can be used as a replacement for Seurat's default normalization methods. To do so, we need to first normalize the data stored in a `SingleCellExperiment` object and then coerce that object to a Seurat object. This can be done with the `as.Seurat` function provided in the `Seurat` package (tested with Seurat 4.0.3).

```{r seurat, eval=FALSE}
library(Seurat)
sce <- PsiNorm(sce)
sce <- logNormCounts(sce)
seu <- as.Seurat(sce)
```

From this point on, one can continue the analysis with the recommended Seurat workflow, but using PsiNorm log-normalized data.

# Using PsiNorm with HDF5 files

Thanks to the `HDF5Array` and `DelayedArray` packages, PsiNorm can be applied directly to HDF5-backed matrices without the need for the user to change the code. As an example, we use a dataset from the `TENxPBMCData` package, which provides several SingleCellExperiment objects with HDF5-backed matrices as their assays.

```{r pbmc, message=FALSE, warning=FALSE}
library(TENxPBMCData)

sce <- TENxPBMCData("pbmc4k")
sce
```

In particular, we use the `pbmc4k` dataset that contains about 4,000 PBMCs from a healthy donor.

The `counts` assay of this object is a `DelayedMatrix` backed by a HDF5 file. Hence, the data are store on disk (out of memory).

```{r seed}
counts(sce)
seed(counts(sce))
```

Thanks to the `DelayedArray` framework, we can apply PsiNorm using the same code that we have used in the case of in-memory data.

```{r psinorm-hdf5}
sce<-PsiNorm(sce)
sce<-logNormCounts(sce)
sce
```

Note that `logNormCounts` is a delayed operation, meaning that the actual log-normalized values will be computed only when needed by the user. In other words, the data are still stored out-of-memory as the original count matrix and the log-normalized data will be computed only when `logcounts(sce)` is realized into memory.

```{r check-seed}
seed(logcounts(sce))
```

# Session Information

```{r}
sessionInfo()
```
