Normalization by distributional resampling of high throughput single-cell RNA-sequencing data
05/30/2021
Package
Dino 1.0.0
Contents
1 Introduction
Over the past decade, advances in single-cell RNA-sequencing (scRNA-seq) technologies have significantly increased the sensitivity and specificity with which cellular transcriptional dynamics can be analyzed. Further, parallel increases in the number cells which can be simultaneously sequenced have allowed for novel analysis pipelines including the description of transcriptional trajectories and the discovery of rare sub-populations of cells. The development of droplet-based, unique-molecular-identifier (UMI) protocols such as Drop-seq, inDrop, and the 10x Genomics Chromium platform have significantly contributed to these advances. In particular, the commercially available 10x Genomics platform has allowed the rapid and cost effective gene expression profiling of hundreds to tens of thousands of cells across many studies to date.
The use of UMIs in the 10x Genomics and related platforms has augmented these developments in sequencing technology by tagging individual mRNA transcripts with unique cell and transcript specific identifiers. In this way, biases due to transcript length and PCR amplification have been significantly reduced. However, technical variability in sequencing depth remains and, consequently, normalization to adjust for sequencing depth is required to ensure accurate downstream analyses. To address this, we introduce Dino and Dino, its corresponding R package.
Dino utilizes a flexible mixture of Negative Binomials model of gene expression to reconstruct full gene-specific expression distributions which are independent of sequencing depth. By giving exact zeros positive probability, the Negative Binomial components are applicable to shallow sequencing (high proportions of zeros). Additionally, the mixture component is robust to cell heterogeneity as it accommodates multiple centers of gene expression in the distribution. By directly modeling (possibly heterogenous) gene-specific expression distributions, Dino outperforms competing approaches, especially for datasets in which the proportion of zeros is high as is typical for modern, UMI based protocols.
Dino does not attempt to correct for batch or other sample specific effects, and will only do so to the extent that they are correlated with sequencing depth. In situations where batch effects are expected, downstream analysis may benefit from such accommodations.
2 Quick Start
2.1 Installation
To install Dino from Github, run
devtools::install_github('JBrownBiostat/Dino')2.2 All-in-one function
Dino (function) is an all-in-one function to normalize raw UMI count data from 10X Cell Ranger or similar protocols. Under default options, Dino outputs a sparse matrix of normalized expression. SeuratFromDino provides one-line functionality to return a Seurat object from raw UMI counts or from a previously normalized expression matrix.
library(Dino)
# Return a sparse matrix of normalized expression
Norm_Mat <- Dino(UMI_Mat)
# Return a Seurat object from already normalized expression
# Use un-transformed normalized expression
Norm_Seurat <- SeuratFromDino(Norm_Mat, doNorm = FALSE, doLog = FALSE)
# Return a Seurat object from UMI expression
# Transform normalized expression as log(x + 1) to improve
# some types of downstream analysis
Norm_Seurat <- SeuratFromDino(UMI_Mat)3 Detailed steps
3.1 Read UMI data
To facilitate concrete examples, we demonstrate normalization on a small subset of sequencing data from about 3,000 peripheral blood mononuclear cells (PBMCs) published by 10X Genomics. This dataset, named pbmcSmall contains 200 cells and 1,000 genes and is included with the Dino package.
set.seed(1)
# Bring pbmcSmall into R environment
library(Dino)
library(Seurat)
library(Matrix)
data("pbmcSmall")
print(dim(pbmcSmall))## [1] 1000 200While Dino was developed to normalize UMI count data, it will run on any matrix of non-negative expression data; user caution is advised if applying Dino to non-UMI sequencing protocols. Input formats may be sparse or dense matrices of expression with genes (features) on the rows and cells (samples) on the columns.
3.2 Clean UMI data
While Dino can normalize the pbmcSmall dataset as it currently exists, the resulting normalized matrix, and in particular, downstream analysis are likely to be improved by cleaning the data. Of greatest use is removing genes that are expected not to contain useful information. This set of genes may be case dependent, but a good rule of thumb for UMI protocols is to remove genes lacking a minimum of non-zero expression prior to normalization and analysis.
By default, Dino will not perform the resampling algorithm on any genes without at least 10 non-zero samples, and will rather normalize such genes by scaling with sequencing depth. To demonstrate a stricter threshold, we remove genes lacking at least 20 non-zero samples prior to normalization.
# Filter genes for a minimum of non-zero expression
pbmcSmall <- pbmcSmall[rowSums(pbmcSmall != 0) >= 20, ]
print(dim(pbmcSmall))## [1] 907 2003.3 Normalize UMI data
Dino contains several options to tune output. One of particular interest is nCores which allows for parallel computation of normalized expression. By default, Dino runs with two threads. Choosing nCores = 0 will utilize all available cores, and otherwise an integer number of parallel instances can be chosen.
# Normalize data
pbmcSmall_Norm <- Dino(pbmcSmall)3.4 Clustering with Seurat
After normalization, Dino makes it easy to perform data analysis. The default output is the normalized matrix in sparse format, and Dino additionally provides a function to transform normalized output into a Seurat object. We demonstrate this by running a quick clustering pipeline in Seurat. Much of the pipeline is modified from the tutorial at https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html
# Reformat normalized expression as a Seurat object
pbmcSmall_Seurat <- SeuratFromDino(pbmcSmall_Norm, doNorm = FALSE)
# Cluster pbmcSmall_Seurat
pbmcSmall_Seurat <- FindVariableFeatures(pbmcSmall_Seurat,
selection.method = "mvp")
pbmcSmall_Seurat <- ScaleData(pbmcSmall_Seurat,
features = rownames(pbmcSmall_Norm))
pbmcSmall_Seurat <- RunPCA(pbmcSmall_Seurat,
features = VariableFeatures(object = pbmcSmall_Seurat),
verbose = FALSE)
pbmcSmall_Seurat <- FindNeighbors(pbmcSmall_Seurat, dims = 1:10)
pbmcSmall_Seurat <- FindClusters(pbmcSmall_Seurat, verbose = FALSE)
pbmcSmall_Seurat <- RunUMAP(pbmcSmall_Seurat, dims = 1:10)
DimPlot(pbmcSmall_Seurat, reduction = "umap")
3.5 Normalizing data formatted as SingleCellExperiment
Dino additionally supports the normalization of datasets formatted as SingleCellExperiment. As with the Seurat pipeline, this functionality is implemented through the use of a wrapper function. We demonstrate this by quickly converting the pbmcSmall dataset to a SingleCellExperiment object and then normalizing.
# Reformatting pbmcSmall as a SingleCellExperiment
library(SingleCellExperiment)
pbmc_SCE <- SingleCellExperiment(assays = list("counts" = pbmcSmall))
# Run Dino
pbmc_SCE <- Dino_SCE(pbmc_SCE)
str(normcounts(pbmc_SCE))## Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## ..@ i : int [1:162768] 0 1 2 3 4 5 6 7 8 9 ...
## ..@ p : int [1:201] 0 812 1617 2423 3233 4041 4856 5672 6468 7271 ...
## ..@ Dim : int [1:2] 907 200
## ..@ Dimnames:List of 2
## .. ..$ : chr [1:907] "ENSG00000087086" "ENSG00000167996" "ENSG00000251562" "ENSG00000205542" ...
## .. ..$ : chr [1:200] "CCAACCTGACGTAC-1" "ATCTGGGATTCCGC-1" "TACTTTCTTTTGGG-1" "CAGGCCGAACACGT-1" ...
## ..@ x : num [1:162768] 104.4 27.8 13.7 30 10.9 ...
## ..@ factors : list()3.6 Alternate sequencing depth
By default, Dino computes sequencing depth, which is corrected for in the normalized data, as the sum of expression for a cell (sample) across genes. This sum is then scaled such that the median depth is 1. For some datasets, however, it may be beneficial to run Dino on an alternately computed set of sequencing depths. Note: it is generally recommended that the median depth not be far from 1 as this corresponds to recomputing expression as though all cells had been sequenced at the median depth.
A simple pipeline to compute alternate sequencing depths utilizes the Scran method for computing normalization scale factors, and is demonstrated below.
library(scran)
# Compute scran size factors
scranSizes <- calculateSumFactors(pbmcSmall)
# Re-normalize data
pbmcSmall_SNorm <- Dino(pbmcSmall, nCores = 1, depth = log(scranSizes))4 Method
4.1 Model
Dino models observed UMI counts as a mixture of Negative Binomial random variables. The Negative Binomial distribution can, however, be decomposed into a hierarchical Gamma-Poisson distribution, so for gene \(g\) and cell \(j\), the Dino model for UMI counts is: \[y_{gj}\sim f^{P}(\lambda_{gj}\delta_{j})\\ \lambda_{gj}\sim\sum_{K}\pi_{k}f^{G}\left(\frac{\mu_{gk}}{\theta_g},\theta_g\right)\] where \(f^{P}\) is a Poisson distribution parameterized by mean \(\lambda_{gj}\delta_{j}\) and \(f^{G}\) is a Gamma distribution parameterized by shape \(\mu_{gk}/\theta_g\) and scale \(\theta_g\). \(\delta_{j}\) is the cell-specific sequencing depth, \(\lambda_{gj}\) is the latent level of gene/cell-specific expression independent of depth, component probabilities \(\pi_{k}\) sum to 1, the Gamma distribution is parameterized such that \(\mu_{gk}\) denotes the distribution mean, and the Gamma scale paramter, \(\theta_g\), is constant across mixture components.
Following model fitting for a fixed gene through an accelerated EM algorithm, Dino produces normalized expression values by resampling from the posterior distribution of the latent expression parameters, \(\lambda_{gj}\). It can be shown that the distribution on the \(\lambda_{j}\) (dropping the gene-specific subscript \(g\) as calculations are repreated across genes) is a mixture of Gammas, specifically: \[\mathbb{P}(\lambda_{j}|y_{j},\delta_j)=\sum_{K}\tau_{kj}f^{G}\left(\frac{\mu_{k}}{\theta}+\gamma y_{j},\frac{1}{\frac{1}{\theta}+\gamma\delta_j}\right)\] where \(\tau_{kj}\) denotes the conditional probability that \(\lambda_{gj}\) was sampled from mixture component \(k\) and \(\gamma\) is a global concentration parameter. The \(\tau_{kj}\) are estimated as part of the implementation of the EM algorithm in Dino. The adjustment from the concentration parameter can be seen as a bias in the normalized values towards a scale-factor version of normalization, since, in the limit of \(\gamma\), the normalized expression for cell \(j\) converges to \(y_j/\delta_j\). Default values of \(\gamma=15\) have proven successful.
4.2 Mixture components \(K\)
Approximating the flexibility of a non-parametric method, Dino uses a large number of mixture components, \(K\), in order to capture the full heterogeneity of expression that may exist for a given gene. The gene-specific number of components is estimated as the square root of the number of strictly positive UMI counts for a given gene. By default, \(K\) is limited to be no larger than 100. In simulation, large values of \(K\) are shown to successfully reconstruct both unimodal and multimodal underlying distributions. For example, when UMI counts are estimated under a single negative binomial distribution, the Dino fitted prior distribution (black, right panel) which is used to sample normalized expression closely matches the theoretical sampling distribution (red, right panel). Likewise, the fitted means (\(\mu_k\) in the model, gray lines, left panel) span the range of the simulated data (heat map of counts, left panel), but concentrate around the theoretical mean of the sampling distribution (red, left panel).
## TableGrob (2 x 2) "arrange": 3 grobs
## z cells name grob
## 1 1 (2-2,1-1) arrange gtable[layout]
## 2 2 (2-2,2-2) arrange gtable[layout]
## 3 3 (1-1,1-2) arrange text[GRID.text.292]Simulating data from a pair of Negative Binomial distributions with different means and different dispersion parameters yields similar results in the multimodal case.
## TableGrob (2 x 2) "arrange": 3 grobs
## z cells name grob
## 1 1 (2-2,1-1) arrange gtable[layout]
## 2 2 (2-2,2-2) arrange gtable[layout]
## 3 3 (1-1,1-2) arrange text[GRID.text.497]5 Session Information
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=C
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggpubr_0.4.0 gridExtra_2.3
## [3] ggplot2_3.3.5 SingleCellExperiment_1.16.0
## [5] SummarizedExperiment_1.24.0 Biobase_2.54.0
## [7] GenomicRanges_1.46.0 GenomeInfoDb_1.30.0
## [9] IRanges_2.28.0 S4Vectors_0.32.0
## [11] BiocGenerics_0.40.0 MatrixGenerics_1.6.0
## [13] matrixStats_0.61.0 Matrix_1.3-4
## [15] SeuratObject_4.0.2 Seurat_4.0.5
## [17] Dino_1.0.0 knitr_1.36
## [19] BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.22
## [3] tidyselect_1.1.1 htmlwidgets_1.5.4
## [5] BiocParallel_1.28.0 Rtsne_0.15
## [7] munsell_0.5.0 ScaledMatrix_1.2.0
## [9] codetools_0.2-18 ica_1.0-2
## [11] statmod_1.4.36 scran_1.22.0
## [13] future_1.22.1 miniUI_0.1.1.1
## [15] withr_2.4.2 colorspace_2.0-2
## [17] highr_0.9 ROCR_1.0-11
## [19] ggsignif_0.6.3 tensor_1.5
## [21] listenv_0.8.0 labeling_0.4.2
## [23] GenomeInfoDbData_1.2.7 polyclip_1.10-0
## [25] farver_2.1.0 parallelly_1.28.1
## [27] vctrs_0.3.8 generics_0.1.1
## [29] xfun_0.27 R6_2.5.1
## [31] rsvd_1.0.5 locfit_1.5-9.4
## [33] bitops_1.0-7 spatstat.utils_2.2-0
## [35] DelayedArray_0.20.0 assertthat_0.2.1
## [37] promises_1.2.0.1 scales_1.1.1
## [39] gtable_0.3.0 beachmat_2.10.0
## [41] globals_0.14.0 goftest_1.2-3
## [43] rlang_0.4.12 splines_4.1.1
## [45] rstatix_0.7.0 lazyeval_0.2.2
## [47] hexbin_1.28.2 spatstat.geom_2.3-0
## [49] broom_0.7.9 BiocManager_1.30.16
## [51] yaml_2.2.1 reshape2_1.4.4
## [53] abind_1.4-5 backports_1.2.1
## [55] httpuv_1.6.3 tools_4.1.1
## [57] bookdown_0.24 ellipsis_0.3.2
## [59] spatstat.core_2.3-0 jquerylib_0.1.4
## [61] RColorBrewer_1.1-2 ggridges_0.5.3
## [63] Rcpp_1.0.7 plyr_1.8.6
## [65] sparseMatrixStats_1.6.0 zlibbioc_1.40.0
## [67] purrr_0.3.4 RCurl_1.98-1.5
## [69] rpart_4.1-15 deldir_1.0-6
## [71] pbapply_1.5-0 cowplot_1.1.1
## [73] zoo_1.8-9 haven_2.4.3
## [75] ggrepel_0.9.1 cluster_2.1.2
## [77] magrittr_2.0.1 data.table_1.14.2
## [79] RSpectra_0.16-0 magick_2.7.3
## [81] scattermore_0.7 openxlsx_4.2.4
## [83] lmtest_0.9-38 RANN_2.6.1
## [85] fitdistrplus_1.1-6 hms_1.1.1
## [87] patchwork_1.1.1 mime_0.12
## [89] evaluate_0.14 xtable_1.8-4
## [91] rio_0.5.27 readxl_1.3.1
## [93] compiler_4.1.1 tibble_3.1.5
## [95] KernSmooth_2.23-20 crayon_1.4.1
## [97] htmltools_0.5.2 mgcv_1.8-38
## [99] later_1.3.0 snow_0.4-3
## [101] tidyr_1.1.4 DBI_1.1.1
## [103] MASS_7.3-54 car_3.0-11
## [105] parallel_4.1.1 metapod_1.2.0
## [107] igraph_1.2.7 forcats_0.5.1
## [109] pkgconfig_2.0.3 foreign_0.8-81
## [111] plotly_4.10.0 scuttle_1.4.0
## [113] spatstat.sparse_2.0-0 bslib_0.3.1
## [115] dqrng_0.3.0 XVector_0.34.0
## [117] stringr_1.4.0 digest_0.6.28
## [119] sctransform_0.3.2 RcppAnnoy_0.0.19
## [121] spatstat.data_2.1-0 cellranger_1.1.0
## [123] rmarkdown_2.11 leiden_0.3.9
## [125] uwot_0.1.10 edgeR_3.36.0
## [127] DelayedMatrixStats_1.16.0 curl_4.3.2
## [129] shiny_1.7.1 lifecycle_1.0.1
## [131] nlme_3.1-153 jsonlite_1.7.2
## [133] carData_3.0-4 BiocNeighbors_1.12.0
## [135] viridisLite_0.4.0 limma_3.50.0
## [137] fansi_0.5.0 pillar_1.6.4
## [139] lattice_0.20-45 fastmap_1.1.0
## [141] httr_1.4.2 survival_3.2-13
## [143] glue_1.4.2 zip_2.2.0
## [145] png_0.1-7 bluster_1.4.0
## [147] stringi_1.7.5 sass_0.4.0
## [149] BiocSingular_1.10.0 dplyr_1.0.7
## [151] irlba_2.3.3 future.apply_1.8.16 Citation
If you use Dino in your analysis, please cite our paper:
Brown, J., Ni, Z., Mohanty, C., Bacher, R., and Kendziorski, C. (2020). “Normalization by distributional resampling of high throughput single-cell RNA-sequencing data.” bioRxiv. https://doi.org/10.1101/2020.10.28.359901.
Other work referenced in this vignette include:
Satija, R., Farrell, J.A., Gennert, D., Schier, A.F. and Regev, A. (2015). “Spatial reconstruction of single-cell gene expression data.” Nat. Biotechnol., 33, 495–502. https://doi.org/10.1038/nbt.3192
Amezquita, R.A., Lun, A.T.L., Becht, E., Carey, V.J., Carpp, L.N., Geistlinger, L., Marini, F., Rue-Albrecht, K., Risso, D., Soneson, C., et al. (2020). “Orchestrating single-cell analysis with Bioconductor.” Nat. Methods, 17, 137–145. https://doi.org/10.1038/s41592-019-0654-x
Lun, A. T. L., Bach, K. and Marioni, J. C. (2016). “Pooling across cells to normalize single-cell RNA sequencing data with many zero counts.” Genome Biol., 17, 1–14. https://doi.org/10.1186/s13059-016-0947-7
7 Contact
Jared Brown: brownj AT biostat DOT wisc DOT edu
Christina Kendziorski: kendzior AT biostat DOT wisc DOT edu