Dimensionality Reduction

Seurat Toolkit

Reduce high-dimensional gene expression data from individual cells into a lower-dimensional space for visualization. This lab explores PCA, tSNE and UMAP.
Authors

Åsa Björklund

Paulo Czarnewski

Susanne Reinsbach

Roy Francis

Published

06-Feb-2024

Note

Code chunks run R commands unless otherwise specified.

1 Data preparation

First, let’s load all necessary libraries and the QC-filtered dataset from the previous step.

# Activate conda environment to get the correct python path
reticulate::use_condaenv("seurat", conda = "/opt/conda/bin/conda")
reticulate::py_discover_config()
python:         /opt/conda/envs/seurat/bin/python
libpython:      /opt/conda/envs/seurat/lib/libpython3.8.so
pythonhome:     /opt/conda/envs/seurat:/opt/conda/envs/seurat
version:        3.8.18 | packaged by conda-forge | (default, Dec 23 2023, 17:21:28)  [GCC 12.3.0]
numpy:          /opt/conda/envs/seurat/lib/python3.8/site-packages/numpy
numpy_version:  1.24.4

NOTE: Python version was forced by use_python function
suppressPackageStartupMessages({
    library(Seurat)
    library(ggplot2) # plotting
    library(patchwork) # combining figures
    library(scran)
})
# download pre-computed data if missing or long compute
fetch_data <- TRUE

# url for source and intermediate data
path_data <- "https://export.uppmax.uu.se/naiss2023-23-3/workshops/workshop-scrnaseq"
path_file <- "data/covid/results/seurat_covid_qc.rds"
if (!dir.exists(dirname(path_file))) dir.create(dirname(path_file), recursive = TRUE)
if (fetch_data && !file.exists(path_file)) download.file(url = file.path(path_data, "covid/results/seurat_covid_qc.rds"), destfile = path_file)
alldata <- readRDS(path_file)

2 Feature selection

We first need to define which features/genes are important in our dataset to distinguish cell types. For this purpose, we need to find genes that are highly variable across cells, which in turn will also provide a good separation of the cell clusters.

suppressWarnings(suppressMessages(alldata <- FindVariableFeatures(alldata, selection.method = "vst", nfeatures = 2000, verbose = FALSE, assay = "RNA")))
top20 <- head(VariableFeatures(alldata), 20)

LabelPoints(plot = VariableFeaturePlot(alldata), points = top20, repel = TRUE)

3 Z-score transformation

Now that the genes have been selected, we now proceed with PCA. Since each gene has a different expression level, it means that genes with higher expression values will naturally have higher variation that will be captured by PCA. This means that we need to somehow give each gene a similar weight when performing PCA (see below). The common practice is to center and scale each gene before performing PCA. This exact scaling called Z-score normalization is very useful for PCA, clustering and plotting heatmaps. Additionally, we can use regression to remove any unwanted sources of variation from the dataset, such as cell cycle, sequencing depth, percent mitochondria etc. This is achieved by doing a generalized linear regression using these parameters as co-variates in the model. Then the residuals of the model are taken as the regressed data. Although perhaps not in the best way, batch effect regression can also be done here. By default, variables are scaled in the PCA step and is not done separately. But it could be achieved by running the commands below:

alldata <- ScaleData(alldata, vars.to.regress = c("percent_mito", "nFeature_RNA"), assay = "RNA")

4 PCA

Performing PCA has many useful applications and interpretations, which much depends on the data used. In the case of single-cell data, we want to segregate samples based on gene expression patterns in the data.

To run PCA, you can use the function RunPCA().

alldata <- RunPCA(alldata, npcs = 50, verbose = F)

We then plot the first principal components.

wrap_plots(
    DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 1:2),
    DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 3:4),
    DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 5:6),
    ncol = 3
) + plot_layout(guides = "collect")

To identify which genes (Seurat) or metadata parameters (Scater/Scran) contribute the most to each PC, one can retrieve the loading matrix information. Unfortunately, this is not implemented in Scater/Scran, so you will need to compute PCA using logcounts.

VizDimLoadings(alldata, dims = 1:5, reduction = "pca", ncol = 5, balanced = T)

We can also plot the amount of variance explained by each PC.

ElbowPlot(alldata, reduction = "pca", ndims = 50)

Based on this plot, we can see that the top 8 PCs retain a lot of information, while other PCs contain progressively less. However, it is still advisable to use more PCs since they might contain information about rare cell types (such as platelets and DCs in this dataset)

5 tSNE

We will now run BH-tSNE.

alldata <- RunTSNE(
    alldata,
    reduction = "pca", dims = 1:30,
    perplexity = 30,
    max_iter = 1000,
    theta = 0.5,
    eta = 200,
    num_threads = 0
)
# see ?Rtsne and ?RunTSNE for more info

We plot the tSNE scatterplot colored by dataset. We can clearly see the effect of batches present in the dataset.

DimPlot(alldata, reduction = "tsne", group.by = "orig.ident")

6 UMAP

We can now run UMAP for cell embeddings.

alldata <- RunUMAP(
    alldata,
    reduction = "pca",
    dims = 1:30,
    n.components = 2,
    n.neighbors = 30,
    n.epochs = 200,
    min.dist = 0.3,
    learning.rate = 1,
    spread = 1
)
# see ?RunUMAP for more info

A feature of UMAP is that it is not limited by the number of dimensions the data cen be reduced into (unlike tSNE). We can simply reduce the dimentions altering the n.components parameter. So here we will create a UMAP with 10 dimensions.

In Seurat, we can add in additional reductions, by default they are named “pca”, “umap”, “tsne” etc. depending on the function you run. Here we will specify an alternative name for the umap with the reduction.name parameter.

alldata <- RunUMAP(
    alldata,
    reduction.name = "UMAP10_on_PCA",
    reduction = "pca",
    dims = 1:30,
    n.components = 10,
    n.neighbors = 30,
    n.epochs = 200,
    min.dist = 0.3,
    learning.rate = 1,
    spread = 1
)
# see ?RunUMAP for more info

UMAP is plotted colored per dataset. Although less distinct as in the tSNE, we still see quite an effect of the different batches in the data.

wrap_plots(
    DimPlot(alldata, reduction = "umap", group.by = "orig.ident") + ggplot2::ggtitle(label = "UMAP_on_PCA"),
    DimPlot(alldata, reduction = "UMAP10_on_PCA", group.by = "orig.ident", dims = 1:2) + ggplot2::ggtitle(label = "UMAP10_on_PCA"),
    DimPlot(alldata, reduction = "UMAP10_on_PCA", group.by = "orig.ident", dims = 3:4) + ggplot2::ggtitle(label = "UMAP10_on_PCA"),
    ncol = 3
) + plot_layout(guides = "collect")

We can now plot PCA, UMAP and tSNE side by side for comparison. Have a look at the UMAP and tSNE. What similarities/differences do you see? Can you explain the differences based on what you learned during the lecture? Also, we can conclude from the dimensionality reductions that our dataset contains a batch effect that needs to be corrected before proceeding to clustering and differential gene expression analysis.

wrap_plots(
    DimPlot(alldata, reduction = "pca", group.by = "orig.ident"),
    DimPlot(alldata, reduction = "tsne", group.by = "orig.ident"),
    DimPlot(alldata, reduction = "umap", group.by = "orig.ident"),
    ncol = 3
) + plot_layout(guides = "collect")

Discuss

We have now done Variable gene selection, PCA and UMAP with the settings we selected for you. Test a few different ways of selecting variable genes, number of PCs for UMAP and check how it influences your embedding.

7 Z-scores & DR graphs

Although running a second dimensionality reduction (i.e tSNE or UMAP) on PCA would be a standard approach (because it allows higher computation efficiency), the options are actually limitless. Below we will show a couple of other common options such as running directly on the scaled data (z-scores) (which was used for PCA) or on a graph built from scaled data. We will only work with UMAPs, but the same applies for tSNE.

7.1 UMAP from z-scores

To run tSNE or UMAP on the scaled data, one first needs to select the number of variables to use. This is because including dimensions that do contribute to the separation of your cell types will in the end mask those differences. Another reason for it is because running with all genes/features also will take longer or might be computationally unfeasible. Therefore we will use the scaled data of the highly variable genes.

alldata <- RunUMAP(
    alldata,
    reduction.name = "UMAP_on_ScaleData",
    features = alldata@assays$RNA@var.features,
    assay = "RNA",
    n.components = 2,
    n.neighbors = 30,
    n.epochs = 200,
    min.dist = 0.3,
    learning.rate = 1,
    spread = 1
)

7.2 UMAP from graph

To run tSNE or UMAP on the a graph, we first need to build a graph from the data. In fact, both tSNE and UMAP first build a graph from the data using a specified distance matrix and then optimize the embedding. Since a graph is just a matrix containing distances from cell to cell and as such, you can run either UMAP or tSNE using any other distance metric desired. Euclidean and Correlation are usually the most commonly used.

# Build Graph
alldata <- FindNeighbors(alldata,
    reduction = "pca",
    assay = "RNA",
    k.param = 20,
    features = alldata@assays$RNA@var.features
)

alldata <- RunUMAP(alldata,
    reduction.name = "UMAP_on_Graph",
    umap.method = "umap-learn",
    graph = "RNA_snn",
    n.epochs = 200,
    assay = "RNA"
)

We can now plot the UMAP comparing both on PCA vs ScaledSata vs Graph.

p1 <- DimPlot(alldata, reduction = "umap", group.by = "orig.ident") + ggplot2::ggtitle(label = "UMAP_on_PCA")
p2 <- DimPlot(alldata, reduction = "UMAP_on_ScaleData", group.by = "orig.ident") + ggplot2::ggtitle(label = "UMAP_on_ScaleData")
p3 <- DimPlot(alldata, reduction = "UMAP_on_Graph", group.by = "orig.ident") + ggplot2::ggtitle(label = "UMAP_on_Graph")
wrap_plots(p1, p2, p3, ncol = 3) + plot_layout(guides = "collect")

8 Genes of interest

Let’s plot some marker genes for different cell types onto the embedding.

Markers Cell Type
CD3E T cells
CD3E CD4 CD4+ T cells
CD3E CD8A CD8+ T cells
GNLY, NKG7 NK cells
MS4A1 B cells
CD14, LYZ, CST3, MS4A7 CD14+ Monocytes
FCGR3A, LYZ, CST3, MS4A7 FCGR3A+ Monocytes
FCER1A, CST3 DCs 
myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7", "FCGR3A", "CST3", "FCER1A")
FeaturePlot(alldata, reduction = "umap", dims = 1:2, features = myfeatures, ncol = 4, order = T) +
    NoLegend() + NoAxes() + NoGrid()

Discuss

Select some of your dimensionality reductions and plot some of the QC stats that were calculated in the previous lab. Can you see if some of the separation in your data is driven by quality of the cells?

9 Save data

We can finally save the object for use in future steps.

saveRDS(alldata, "data/covid/results/seurat_covid_qc_dr.rds")

10 Session info

Click here 
sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 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: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] scran_1.30.0                scuttle_1.12.0             
 [3] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
 [5] Biobase_2.62.0              GenomicRanges_1.54.1       
 [7] GenomeInfoDb_1.38.5         IRanges_2.36.0             
 [9] S4Vectors_0.40.2            BiocGenerics_0.48.1        
[11] MatrixGenerics_1.14.0       matrixStats_1.0.0          
[13] patchwork_1.1.2             ggplot2_3.4.2              
[15] SeuratObject_4.1.3          Seurat_4.3.0               

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3        rstudioapi_0.14          
  [3] jsonlite_1.8.5            magrittr_2.0.3           
  [5] spatstat.utils_3.0-3      farver_2.1.1             
  [7] rmarkdown_2.22            zlibbioc_1.48.0          
  [9] vctrs_0.6.2               ROCR_1.0-11              
 [11] DelayedMatrixStats_1.24.0 spatstat.explore_3.2-1   
 [13] RCurl_1.98-1.12           S4Arrays_1.2.0           
 [15] htmltools_0.5.5           BiocNeighbors_1.20.2     
 [17] SparseArray_1.2.3         sctransform_0.3.5        
 [19] parallelly_1.36.0         KernSmooth_2.23-20       
 [21] htmlwidgets_1.6.2         ica_1.0-3                
 [23] plyr_1.8.8                plotly_4.10.2            
 [25] zoo_1.8-12                igraph_1.4.3             
 [27] mime_0.12                 lifecycle_1.0.3          
 [29] pkgconfig_2.0.3           rsvd_1.0.5               
 [31] Matrix_1.5-4              R6_2.5.1                 
 [33] fastmap_1.1.1             GenomeInfoDbData_1.2.11  
 [35] fitdistrplus_1.1-11       future_1.32.0            
 [37] shiny_1.7.4               digest_0.6.31            
 [39] colorspace_2.1-0          rprojroot_2.0.3          
 [41] tensor_1.5                dqrng_0.3.0              
 [43] irlba_2.3.5.1             beachmat_2.18.0          
 [45] labeling_0.4.2            progressr_0.13.0         
 [47] fansi_1.0.4               spatstat.sparse_3.0-1    
 [49] httr_1.4.6                polyclip_1.10-4          
 [51] abind_1.4-5               compiler_4.3.0           
 [53] here_1.0.1                withr_2.5.0              
 [55] BiocParallel_1.36.0       MASS_7.3-58.4            
 [57] DelayedArray_0.28.0       bluster_1.12.0           
 [59] tools_4.3.0               lmtest_0.9-40            
 [61] httpuv_1.6.11             future.apply_1.11.0      
 [63] goftest_1.2-3             glue_1.6.2               
 [65] nlme_3.1-162              promises_1.2.0.1         
 [67] grid_4.3.0                Rtsne_0.16               
 [69] cluster_2.1.4             reshape2_1.4.4           
 [71] generics_0.1.3            gtable_0.3.3             
 [73] spatstat.data_3.0-1       tidyr_1.3.0              
 [75] data.table_1.14.8         metapod_1.10.1           
 [77] ScaledMatrix_1.10.0       BiocSingular_1.18.0      
 [79] sp_1.6-1                  utf8_1.2.3               
 [81] XVector_0.42.0            spatstat.geom_3.2-1      
 [83] RcppAnnoy_0.0.20          ggrepel_0.9.3            
 [85] RANN_2.6.1                pillar_1.9.0             
 [87] stringr_1.5.0             limma_3.58.1             
 [89] later_1.3.1               splines_4.3.0            
 [91] dplyr_1.1.2               lattice_0.21-8           
 [93] survival_3.5-5            deldir_1.0-9             
 [95] tidyselect_1.2.0          locfit_1.5-9.8           
 [97] miniUI_0.1.1.1            pbapply_1.7-0            
 [99] knitr_1.43                gridExtra_2.3            
[101] edgeR_4.0.7               scattermore_1.2          
[103] xfun_0.39                 statmod_1.5.0            
[105] stringi_1.7.12            lazyeval_0.2.2           
[107] yaml_2.3.7                evaluate_0.21            
[109] codetools_0.2-19          tibble_3.2.1             
[111] cli_3.6.1                 uwot_0.1.14              
[113] xtable_1.8-4              reticulate_1.30          
[115] munsell_0.5.0             Rcpp_1.0.10              
[117] globals_0.16.2            spatstat.random_3.1-5    
[119] png_0.1-8                 parallel_4.3.0           
[121] ellipsis_0.3.2            sparseMatrixStats_1.14.0 
[123] bitops_1.0-7              listenv_0.9.0            
[125] viridisLite_0.4.2         scales_1.2.1             
[127] ggridges_0.5.4            crayon_1.5.2             
[129] leiden_0.4.3              purrr_1.0.1              
[131] rlang_1.1.1               cowplot_1.1.1