Normalisation
NadineBestard
10/07/2022
Set-up
library(here) # for reproducible paths
library(SingleCellExperiment)
library(scater) # For qc and visualisation
library(scran) # For normalisation
library(Matrix) # For log transorming the raw data
library(ggplot2) # To add titles to plots# Adapted function from VISION to log tranform sparse matrix
# I could not download the package
matLog2 <- function(spmat, scale = FALSE, scaleFactor = 1e6) {
if (scale == TRUE) {
spmat <- t( t(spmat) / colSums(spmat)) * scaleFactor
}
if (is(spmat, "sparseMatrix")) {
matsum <- summary(spmat)
logx <- log2(matsum$x + 1)
logmat <- sparseMatrix(i = matsum$i, j = matsum$j,
x = logx, dims = dim(spmat),
dimnames = dimnames(spmat))
} else {
logmat <- log2(spmat + 1)
}
return(logmat)
}project <- "fire-mice"Normalisation by deconvolution
In order to correct for systematic differences in sequencing coverage between libraries we will normalise the dataset. This involves dividing all counts for each cell by a cell-specific scaling factor, often called a “size factor” (Anders and Huber 2010). The assumption here is that any cell-specific bias (e.g., in capture or amplification efficiency) affects all genes equally via scaling of the expected mean count for that cell. The size factor for each cell represents the estimate of the relative bias in that cell, so division of its counts by its size factor should remove that bias.
Specifically we will used the deconvolution method available in the scran package. This method allows to take in consideration the composition bias between samples (Lun et al., 2016)
# Only compute if first time
if (!(file.exists(here("processed", project, "sce_norm_01.RDS")))) {
sce <- readRDS(here("processed", project, "sce_QC_01.RDS"))
# For reproducibility
set.seed(100)
# Quick clustering to pool samples together and deal with 0 counts
quick_clusters <- quickCluster(sce)
# Calculate size factors
sce <- computeSumFactors(sce, cluster = quick_clusters, min.mean = 0.1)
# Check that there are not negative size factors
summary(sizeFactors(sce))
# Apply size factors and log transform them
sce <- logNormCounts(sce)
# Also log normalise the raw counts
assay(sce, "logcounts_raw") <- matLog2(counts(sce))
saveRDS(sce, here("processed", project, "sce_norm_01.RDS"))
} else{
sce <- readRDS(here("processed", project, "sce_norm_01.RDS"))
}On top of normalisation the data is also log-transformed. The log-transformation is useful as differences in the log-values represent log-fold changes in expression. Or in other words, which is more interesting - a gene that is expressed at an average count of 50 in cell type A and 10 in cell type B, or a gene that is expressed at an average count of 1100 in A and 1000 in B? Log-transformation focuses on the former by promoting contributions from genes with strong relative differences.
Assess Confunding factors impact
Variance Explained plots
Variable-level metrics are computed by the getVarianceExplained() function (before and after normalization). This calculates the percentage of variance of each gene’s expression that is explained by each variable in the colData of the SingleCellExperiment object. We can then use this to determine which experimental factors are contributing most to the variance in expression. This is useful for diagnosing batch effects or to quickly verify that a treatment has an effect.
The percentage of variance explained by a factor is on the x axis, and in the y axis there is the density of the R-squared values across all genes.
The “total” label is the total number of molecules, that correlates with the detected number of genes, “detected”.
Before normalisation
Before normalisation it is expected that most variance will be explained by the sequencing depth, i.e. the total number of umis and the total number of genes
# Before normalisation
# Only compute if first time
if (!(file.exists(here("processed", project, "variance_explained.RDS")))) {
# Calculate the matrix (time consuming step)
var <- getVarianceExplained(
sce,
exprs_values = "logcounts_raw",
variables = c(
"chip",
"genotype",
"mouse",
"subsets_mt_percent",
"detected",
"total"
)
)
saveRDS(var, here("processed", project, "variance_explained.RDS"))
#If not just load created object
} else {
var <- readRDS(here("processed", project, "variance_explained.RDS"))
}
plotExplanatoryVariables(var)
After normalisation
We can see how there is less variance explained now by factors such as the detected genes or the number of counts
# After normalisation
if (!(file.exists(here("processed", project, "variance_explained_norm.RDS")
))) {
var_norm <- getVarianceExplained(
sce,
variables = c(
"chip",
"genotype",
"mouse",
"subsets_mt_percent",
"detected",
"total"
)
)
saveRDS(var_norm, here("processed", project, "variance_explained_norm.RDS"))
} else{
var_norm <- readRDS(here("processed", project, "variance_explained_norm.RDS"))
}
plotExplanatoryVariables(var_norm)## Warning in self$trans$transform(x): NaNs produced## Warning: Transformation introduced infinite values in continuous x-axis## Warning: Removed 1 rows containing non-finite values (stat_density).
Dimensional reduction
We will more accurate dimensional reductions in the next step, only using the most variable genes to reduce noise Another way to assess the variance is with a PCA plot. Here again we can see how the sequencing depth(sum) explains most of the variance before the normalisation
raw <- runPCA(sce, exprs_values = "logcounts_raw")
plotPCA(raw, colour_by= "chip", size_by="sum") + ggtitle("Before normalisation")
sce <- runPCA(sce)
plotPCA(sce, colour_by= "chip", size_by="sum") + ggtitle("After normalisation")
plotPCA(sce, colour_by= "chip", point_size=0.1) +
ggtitle("After normalisation, small dots")
Another type of dimensional reduction are the non linear UMAP and TSNE reductions.
sce <- runUMAP(sce, dimred="PCA")
plotReducedDim(sce, colour_by= "chip", point_size=0.1, dimred = "UMAP") +
ggtitle("UMAP dimensional reduction")
sce <- runTSNE(sce, dimred="PCA")
plotReducedDim(sce, colour_by= "chip", point_size=0.1, dimred = "TSNE") +
ggtitle("TSNE dimensional reduction")
Session Info
Click to expand
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252
## [2] LC_CTYPE=English_United Kingdom.1252
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.1252
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] Matrix_1.4-0 scran_1.22.1
## [3] scater_1.23.5 ggplot2_3.3.5
## [5] scuttle_1.4.0 SingleCellExperiment_1.16.0
## [7] SummarizedExperiment_1.24.0 Biobase_2.54.0
## [9] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1
## [11] IRanges_2.28.0 S4Vectors_0.32.3
## [13] BiocGenerics_0.40.0 MatrixGenerics_1.6.0
## [15] matrixStats_0.61.0 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 RcppAnnoy_0.0.19
## [3] rprojroot_2.0.2 tools_4.1.1
## [5] bslib_0.3.1 utf8_1.2.2
## [7] R6_2.5.1 irlba_2.3.5
## [9] vipor_0.4.5 uwot_0.1.11
## [11] DBI_1.1.2 colorspace_2.0-2
## [13] withr_2.4.3 tidyselect_1.1.1
## [15] gridExtra_2.3 compiler_4.1.1
## [17] cli_3.2.0 BiocNeighbors_1.12.0
## [19] DelayedArray_0.20.0 labeling_0.4.2
## [21] sass_0.4.0 scales_1.1.1
## [23] stringr_1.4.0 digest_0.6.29
## [25] rmarkdown_2.11 XVector_0.34.0
## [27] pkgconfig_2.0.3 htmltools_0.5.2
## [29] sparseMatrixStats_1.6.0 highr_0.9
## [31] fastmap_1.1.0 limma_3.50.0
## [33] rlang_1.0.1 rstudioapi_0.13
## [35] DelayedMatrixStats_1.16.0 farver_2.1.0
## [37] jquerylib_0.1.4 generics_0.1.2
## [39] jsonlite_1.7.3 BiocParallel_1.28.3
## [41] dplyr_1.0.8 RCurl_1.98-1.6
## [43] magrittr_2.0.2 BiocSingular_1.10.0
## [45] GenomeInfoDbData_1.2.7 Rcpp_1.0.8
## [47] ggbeeswarm_0.6.0 munsell_0.5.0
## [49] fansi_1.0.2 viridis_0.6.2
## [51] lifecycle_1.0.1 stringi_1.7.6
## [53] yaml_2.2.2 edgeR_3.36.0
## [55] zlibbioc_1.40.0 Rtsne_0.15
## [57] grid_4.1.1 parallel_4.1.1
## [59] ggrepel_0.9.1 dqrng_0.3.0
## [61] crayon_1.5.0 lattice_0.20-45
## [63] cowplot_1.1.1 beachmat_2.10.0
## [65] locfit_1.5-9.4 metapod_1.2.0
## [67] knitr_1.37 pillar_1.7.0
## [69] igraph_1.2.11 codetools_0.2-18
## [71] ScaledMatrix_1.2.0 glue_1.6.1
## [73] evaluate_0.14 vctrs_0.3.8
## [75] gtable_0.3.0 purrr_0.3.4
## [77] assertthat_0.2.1 xfun_0.29
## [79] rsvd_1.0.5 RSpectra_0.16-0
## [81] viridisLite_0.4.0 tibble_3.1.6
## [83] beeswarm_0.4.0 cluster_2.1.2
## [85] bluster_1.4.0 statmod_1.4.36
## [87] ellipsis_0.3.2