Dimensionality reduction

Material

Exercises

Load the seu dataset you have created yesterday:

seu <- readRDS("day2/seu_day2-1.rds")

And load the following packages:

Once the data is normalized, scaled and variable features have been identified, we can start to reduce the dimensionality of the data. For the PCA, by default, only the previously determined variable features are used as input, but can be defined using features argument if you wish to specify a vector of genes. The PCA will only be run on the variable features, that you can check with VariableFeatures(seu).

seu <- Seurat::RunPCA(seu)

To view the PCA plot:

Seurat::DimPlot(seu, reduction = "pca")

We can colour the PCA plot according to any factor that is present in @meta.data, or for any gene. For example we can take the column percent.globin:

Seurat::FeaturePlot(seu, reduction = "pca", features = "percent.globin")

Note

Note that we used a different plotting function here: FeaturePlot. The difference between DimPlot and FeaturePlot is that the first allows you to color the points in the plot according to a grouping variable (e.g. sample) while the latter allows you to color the points according to a continuous variable (e.g. gene expression).

Exercise

Generate a PCA plot where color is according to counts of a gene (i.e. gene expression). For example, you can take HBA1 (alpha subunit of hemoglobin), or one of the most variable genes (e.g. IGKC).

Generating a PCA plot coloured according to gene expression (here HBA1):

Seurat::FeaturePlot(seu, reduction = "pca", features = "HBA1")

We can generate heatmaps according to their principal component scores calculated in the rotation matrix:

Seurat::DimHeatmap(seu, dims = 1:12, cells = 500, balanced = TRUE)

The elbowplot can help you in determining how many PCs to use for downstream analysis such as UMAP:

Seurat::ElbowPlot(seu, ndims = 40) +
  ggplot2::geom_vline(xintercept = 25, linetype = "dashed", color = "red", linewidth = 1)

The elbow plot ranks principle components based on the percentage of variance explained by each one. Where we observe an “elbow” or flattening curve, the majority of true signal is captured by this number of PCs, eg around 25 PCs for the seu dataset.

Including too many PCs usually does not affect much the result, while including too few PCs can affect the results very much.

UMAP: The goal of these algorithms is to learn the underlying manifold of the data in order to place similar cells together in low-dimensional space.

seu <- Seurat::RunUMAP(seu, dims = 1:25)

To view the UMAP plot:

Seurat::DimPlot(seu, reduction = "umap")

Exercise

Try to change:

A. Color the dots in the UMAP according to a variable (e.g. percent.globin or HBA1). Any idea where the erythrocytes probably are in the UMAP?

B. The number of neighbors used for the calculation of the UMAP. Which is the parameter to change and how did it affect the output. What is the default ? In which situation would you lower/increase this ?

C. The number of dims to extremes dims = 1:5 or dims = 1:50 how did it affect the output ? In your opinion better few PCAs too much or too few ? Why does dims = 1:100 not work? When would more precision be needed?

3D UMAP

Viewing data in 3D could reveal hidden clusters, trajectories, or sample-specific patterns that are not visible in 2D UMAPs — especially in complex single-cell datasets.

library(plotly)
seu <- Seurat::RunUMAP(seu, dims = 1:25, n.components = 3, reduction.name = "umap_3D")

plot_3d_umap <- function(seu,
                         reduction = "umap_3D",
                         sample_col = "orig.ident",
                         point_size = 4,
                         point_opacity = 0.5) {
  
  ## 1. Extract coordinates -------------------------------------------------
  umap_coords <- as.data.frame(Embeddings(seu, reduction = reduction))
  umap_coords$Sample <- seu@meta.data[[sample_col]]
  umap_coords$Cell   <- rownames(umap_coords)   # for hover
  
  ## 2. Build Plotly object --------------------------------------------------
  p <- plot_ly(
    data = umap_coords,
    x = ~umap3D_1, y = ~umap3D_2, z = ~umap3D_3,
    color = ~Sample,                     # default Plotly palette
    type = "scatter3d",
    mode = "markers",
    marker = list(
      size = point_size,
      opacity = point_opacity
    ),
    text = ~paste("Cell:", Cell, "<br>Sample:", Sample),
    hovertemplate = "<b>%{text}</b><extra></extra>"
  ) %>%
    layout(
      title = list(
        text = "<b>3D UMAP of Cells by Sample</b>",
        font = list(size = 18, family = "Arial, sans-serif"),
        x = 0.5, xanchor = "center", y = 0.95
      ),
      scene = list(
        xaxis = list(
          title = "<b>UMAP 1</b>", titlefont = list(size = 14),
          showgrid = TRUE, gridcolor = "rgba(200,200,200,0.5)",
          backgroundcolor = "rgba(245,245,255,0.9)",
          zerolinecolor = "rgba(100,100,100,0.3)"
        ),
        yaxis = list(
          title = "<b>UMAP 2</b>", titlefont = list(size = 14),
          showgrid = TRUE, gridcolor = "rgba(200,200,200,0.5)",
          backgroundcolor = "rgba(245,245,255,0.9)",
          zerolinecolor = "rgba(100,100,100,0.3)"
        ),
        zaxis = list(
          title = "<b>UMAP 3</b>", titlefont = list(size = 14),
          showgrid = TRUE, gridcolor = "rgba(200,200,200,0.5)",
          backgroundcolor = "rgba(245,245,255,0.9)",
          zerolinecolor = "rgba(100,100,100,0.3)"
        ),
        camera = list(eye = list(x = 1.5, y = 1.5, z = 1.5)),
        aspectmode = "cube"
      ),
      legend = list(
        title = list(text = "<b>Sample</b>", font = list(size = 12)),
        bgcolor = "rgba(255,255,255,0.8)",
        bordercolor = "gray", borderwidth = 1
      ),
      margin = list(l = 0, r = 0, t = 60, b = 0),
      paper_bgcolor = "rgba(250,250,252,1)",
      plot_bgcolor  = "rgba(250,250,252,1)"
    ) %>%
    config(
      displayModeBar = TRUE,
      displaylogo = FALSE,
      modeBarButtonsToRemove = c("sendDataToCloud", "lasso2d", "select2d")
    )
  
  return(p)
}

plot_3d_umap(seu = seu, reduction = "umap_3D", sample_col = "orig.ident")

Answer A

Seurat::FeaturePlot(seu, features = c("HBA1",
                                      "percent.globin",
                                      "IGKC",
                                      "percent.mito"))

The erythrocytes are probably in the cluster with a higher percentage of globin expression.

Answer B

seu <- Seurat::RunUMAP(seu, dims = 1:25, n.neighbors = 5)
Seurat::DimPlot(seu, reduction = "umap")

Also in 3D

seu <- Seurat::RunUMAP(seu, dims = 1:25, n.components = 3, n.neighbors = 5, reduction.name = "umap_3D")
13:55:07 UMAP embedding parameters a = 0.9922 b = 1.112
13:55:07 Read 6830 rows and found 25 numeric columns
13:55:07 Using Annoy for neighbor search, n_neighbors = 5
13:55:07 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:55:07 Writing NN index file to temp file /var/folders/cm/nwh7kr5n65b3m40y__9_t5hnkkqk0y/T//RtmpuQDqyu/file5db755e1c849
13:55:07 Searching Annoy index using 1 thread, search_k = 500
13:55:08 Annoy recall = 100%
13:55:08 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 5
13:55:08 Found 3 connected components, falling back to 'spca' initialization with init_sdev = 1
13:55:08 Using 'irlba' for PCA
13:55:08 PCA: 3 components explained 57.33% variance
13:55:08 Scaling init to sdev = 1
13:55:08 Commencing optimization for 500 epochs, with 42540 positive edges
13:55:08 Using rng type: pcg
13:55:11 Optimization finished
plot_3d_umap(seu = seu, reduction = "umap_3D", sample_col = "orig.ident")

The default number of neighbours is 30. It can be of interest to change the number of neighbors if one has subset the data (for instance in the situation where you would only consider the t-cells in your data set), then maybe the number of neighbors in a cluster would anyway be most of the time lower than 30 then 30 is too much. In the other extreme where your dataset is extremely big an increase in the number of neighbors can be considered.

Answer C

seu <- Seurat::RunUMAP(seu, dims = 1:5)
Seurat::DimPlot(seu, reduction = "umap")

seu <- Seurat::RunUMAP(seu, dims = 1:50) 
Seurat::DimPlot(seu, reduction = "umap")

Taking dims = 1:100 does not work as in the step RunPCA by default only 50pcs are calculated, so the maximum that we can consider in further steps are 50, if more precision makes sense, for instance, if the genes that is of interest for your study is not present when the RunPCA was calculated, then an increase in the number of components calculated at start might be interesting tobe considered. Taking too few PCs we have a « blob » everything looks connected. Too many PCs tends to separate everything. Personally it is more interesting for me too have maybe 2 clusters separated of epithelial cells that I then group for further downstream analysis rather than having very distinct cells being clustered together. So I would rather take the « elbow » of the elbow plot a bit further to the right.

After having done these exercises, change the UMAP back to a UMAP based on the first 25 PCs, in order to replicate the exercises in the following chapters. Do this by:

seu <- Seurat::RunUMAP(seu, dims = 1:25)
seu <- Seurat::RunUMAP(seu, dims = 1:25, n.components = 3, reduction.name = "umap_3D")
saveRDS(seu, "day2/seu_day2-2.rds")