Contents

1 Introduction

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.

2 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.

3 Data Simulation

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

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.

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

4 PsiNorm data normalization

5 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.

sce<-PsiNorm(sce)
sce<-logNormCounts(sce)
head(sizeFactors(sce))
#>     Cell1     Cell2     Cell3     Cell4     Cell5     Cell6 
#> 1.1017454 0.9667386 1.0169251 0.9343382 1.0966308 1.1845135

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.

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.

5.1 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.

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)
ARI from raw counts: 0.295
ARI from PsiNorm normalized data: 0.788

Pareto normalization considerably increases the ARI index.

6 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.

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)
Silhouette from raw counts: 0.168
Silhouette from PsiNorm normalized data: 0.534

Pareto normalization considerably increases the Silhouette index.

7 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.

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)
The Correlation with the raw data is: 0.269
The Correlation with the PsiNorm normalized data is: 0.189

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

8 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.

9 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).

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.

10 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.

library(TENxPBMCData)

sce <- TENxPBMCData("pbmc4k")
sce
#> class: SingleCellExperiment 
#> dim: 33694 4340 
#> metadata(0):
#> assays(1): counts
#> rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
#>   ENSG00000268674
#> rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
#> colnames: NULL
#> colData names(11): Sample Barcode ... Individual Date_published
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):

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).

counts(sce)
#> <33694 x 4340> sparse DelayedMatrix object of type "integer":
#>                    [,1]    [,2]    [,3]    [,4] ... [,4337] [,4338] [,4339]
#> ENSG00000243485       0       0       0       0   .       0       0       0
#> ENSG00000237613       0       0       0       0   .       0       0       0
#> ENSG00000186092       0       0       0       0   .       0       0       0
#> ENSG00000238009       0       0       0       0   .       0       0       0
#> ENSG00000239945       0       0       0       0   .       0       0       0
#>             ...       .       .       .       .   .       .       .       .
#> ENSG00000277856       0       0       0       0   .       0       0       0
#> ENSG00000275063       0       0       0       0   .       0       0       0
#> ENSG00000271254       0       0       0       0   .       0       0       0
#> ENSG00000277475       0       0       0       0   .       0       0       0
#> ENSG00000268674       0       0       0       0   .       0       0       0
#>                 [,4340]
#> ENSG00000243485       0
#> ENSG00000237613       0
#> ENSG00000186092       0
#> ENSG00000238009       0
#> ENSG00000239945       0
#>             ...       .
#> ENSG00000277856       0
#> ENSG00000275063       0
#> ENSG00000271254       0
#> ENSG00000277475       0
#> ENSG00000268674       0
seed(counts(sce))
#> An object of class "HDF5ArraySeed"
#> Slot "filepath":
#> [1] "/home/biocbuild/.cache/R/ExperimentHub/bb5702a4107d7_1611"
#> 
#> Slot "name":
#> [1] "/counts"
#> 
#> Slot "as_sparse":
#> [1] TRUE
#> 
#> Slot "type":
#> [1] NA
#> 
#> Slot "dim":
#> [1] 33694  4340
#> 
#> Slot "chunkdim":
#> [1] 512  66
#> 
#> Slot "first_val":
#> [1] 0

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

sce<-PsiNorm(sce)
sce<-logNormCounts(sce)
sce
#> class: SingleCellExperiment 
#> dim: 33694 4340 
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
#>   ENSG00000268674
#> rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
#> colnames: NULL
#> colData names(12): Sample Barcode ... Date_published sizeFactor
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):

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.

seed(logcounts(sce))
#> An object of class "HDF5ArraySeed"
#> Slot "filepath":
#> [1] "/home/biocbuild/.cache/R/ExperimentHub/bb5702a4107d7_1611"
#> 
#> Slot "name":
#> [1] "/counts"
#> 
#> Slot "as_sparse":
#> [1] TRUE
#> 
#> Slot "type":
#> [1] NA
#> 
#> Slot "dim":
#> [1] 33694  4340
#> 
#> Slot "chunkdim":
#> [1] 512  66
#> 
#> Slot "first_val":
#> [1] 0

11 Session Information

sessionInfo()
#> R version 4.4.0 beta (2024-04-15 r86425)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] TENxPBMCData_1.21.0         HDF5Array_1.32.0           
#>  [3] rhdf5_2.48.0                DelayedArray_0.30.0        
#>  [5] SparseArray_1.4.0           S4Arrays_1.4.0             
#>  [7] abind_1.4-5                 Matrix_1.7-0               
#>  [9] scone_1.28.0                cluster_2.1.6              
#> [11] scater_1.32.0               ggplot2_3.5.1              
#> [13] scuttle_1.14.0              splatter_1.28.0            
#> [15] SingleCellExperiment_1.26.0 SummarizedExperiment_1.34.0
#> [17] Biobase_2.64.0              GenomicRanges_1.56.0       
#> [19] GenomeInfoDb_1.40.0         IRanges_2.38.0             
#> [21] S4Vectors_0.42.0            BiocGenerics_0.50.0        
#> [23] MatrixGenerics_1.16.0       matrixStats_1.3.0          
#> [25] BiocStyle_2.32.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] splines_4.4.0             BiocIO_1.14.0            
#>   [3] bitops_1.0-7              filelock_1.0.3           
#>   [5] tibble_3.2.1              R.oo_1.26.0              
#>   [7] XML_3.99-0.16.1           lifecycle_1.0.4          
#>   [9] httr2_1.0.1               pwalign_1.0.0            
#>  [11] edgeR_4.2.0               prabclus_2.3-3           
#>  [13] lattice_0.22-6            MASS_7.3-60.2            
#>  [15] backports_1.4.1           magrittr_2.0.3           
#>  [17] plotly_4.10.4             limma_3.60.0             
#>  [19] sass_0.4.9                rmarkdown_2.26           
#>  [21] jquerylib_0.1.4           yaml_2.3.8               
#>  [23] flexmix_2.3-19            cowplot_1.1.3            
#>  [25] bayesm_3.1-6              DBI_1.2.2                
#>  [27] RColorBrewer_1.1-3        ShortRead_1.62.0         
#>  [29] zlibbioc_1.50.0           purrr_1.0.2              
#>  [31] mixtools_2.0.0            R.utils_2.12.3           
#>  [33] RCurl_1.98-1.14           nnet_7.3-19              
#>  [35] tensorA_0.36.2.1          rappdirs_0.3.3           
#>  [37] GenomeInfoDbData_1.2.12   ggrepel_0.9.5            
#>  [39] irlba_2.3.5.1             RSpectra_0.16-1          
#>  [41] svglite_2.1.3             DelayedMatrixStats_1.26.0
#>  [43] codetools_0.2-20          xml2_1.3.6               
#>  [45] tidyselect_1.2.1          farver_2.1.1             
#>  [47] UCSC.utils_1.0.0          ScaledMatrix_1.12.0      
#>  [49] viridis_0.6.5             BiocFileCache_2.12.0     
#>  [51] GenomicAlignments_1.40.0  jsonlite_1.8.8           
#>  [53] BiocNeighbors_1.22.0      survival_3.6-4           
#>  [55] systemfonts_1.0.6         segmented_2.0-4          
#>  [57] tools_4.4.0               progress_1.2.3           
#>  [59] rARPACK_0.11-0            Rcpp_1.0.12              
#>  [61] glue_1.7.0                gridExtra_2.3            
#>  [63] xfun_0.43                 dplyr_1.1.4              
#>  [65] withr_3.0.0               BiocManager_1.30.22      
#>  [67] fastmap_1.1.1             rhdf5filters_1.16.0      
#>  [69] latticeExtra_0.6-30       boot_1.3-30              
#>  [71] fansi_1.0.6               caTools_1.18.2           
#>  [73] digest_0.6.35             rsvd_1.0.5               
#>  [75] mime_0.12                 R6_2.5.1                 
#>  [77] colorspace_2.1-0          gtools_3.9.5             
#>  [79] jpeg_0.1-10               biomaRt_2.60.0           
#>  [81] RSQLite_2.3.6             diptest_0.77-1           
#>  [83] R.methodsS3_1.8.2         tidyr_1.3.1              
#>  [85] hexbin_1.28.3             utf8_1.2.4               
#>  [87] generics_0.1.3            data.table_1.15.4        
#>  [89] rtracklayer_1.64.0        robustbase_0.99-2        
#>  [91] class_7.3-22              htmlwidgets_1.6.4        
#>  [93] prettyunits_1.2.0         httr_1.4.7               
#>  [95] pkgconfig_2.0.3           gtable_0.3.5             
#>  [97] modeltools_0.2-23         blob_1.2.4               
#>  [99] hwriter_1.3.2.1           XVector_0.44.0           
#> [101] htmltools_0.5.8.1         bookdown_0.39            
#> [103] kableExtra_1.4.0          scales_1.3.0             
#> [105] png_0.1-8                 rstudioapi_0.16.0        
#> [107] knitr_1.46                rjson_0.2.21             
#> [109] nlme_3.1-164              checkmate_2.3.1          
#> [111] curl_5.2.1                cachem_1.0.8             
#> [113] stringr_1.5.1             BiocVersion_3.19.1       
#> [115] KernSmooth_2.23-22        parallel_4.4.0           
#> [117] vipor_0.4.7               RUVSeq_1.38.0            
#> [119] AnnotationDbi_1.66.0      restfulr_0.0.15          
#> [121] pillar_1.9.0              grid_4.4.0               
#> [123] vctrs_0.6.5               gplots_3.1.3.1           
#> [125] BiocSingular_1.20.0       dbplyr_2.5.0             
#> [127] beachmat_2.20.0           beeswarm_0.4.0           
#> [129] evaluate_0.23             magick_2.8.3             
#> [131] tinytex_0.50              GenomicFeatures_1.56.0   
#> [133] cli_3.6.2                 locfit_1.5-9.9           
#> [135] compiler_4.4.0            Rsamtools_2.20.0         
#> [137] rlang_1.1.3               crayon_1.5.2             
#> [139] labeling_0.4.3            mclust_6.1.1             
#> [141] interp_1.1-6              aroma.light_3.34.0       
#> [143] ggbeeswarm_0.7.2          stringi_1.8.3            
#> [145] viridisLite_0.4.2         deldir_2.0-4             
#> [147] BiocParallel_1.38.0       munsell_0.5.1            
#> [149] Biostrings_2.72.0         lazyeval_0.2.2           
#> [151] EDASeq_2.38.0             compositions_2.0-8       
#> [153] ExperimentHub_2.12.0      hms_1.1.3                
#> [155] sparseMatrixStats_1.16.0  bit64_4.0.5              
#> [157] Rhdf5lib_1.26.0           KEGGREST_1.44.0          
#> [159] fpc_2.2-12                statmod_1.5.0            
#> [161] AnnotationHub_3.12.0      highr_0.10               
#> [163] kernlab_0.9-32            memoise_2.0.1            
#> [165] bslib_0.7.0               DEoptimR_1.1-3           
#> [167] bit_4.0.5