Preparation

Install packages

You need to do this only once!
The script below checks whether the listed CRAN packages are installed on your computer. If they are not, they will be automatically installed.

## First specify the packages of interest
packages = c("devtools","hdf5r","dplyr","ggplot2","stringr",
             "RColorBrewer","useful","readr","BiocManager")

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)
install.packages("Matrix")

Seurat for scRNAseq data analysis

install.packages('Seurat')

The following is necessary for interaction with hdf5 objects, and can only be installed using BiocManager installer.

# You need this only once
BiocManager::install("rhdf5")

1. Introduction

In this COO, you will learn how to analyse single cell RNA sequencing data (scRNAseq) using the well established Seurat pipeline.

This pipeline includes functions that do most of the work ‘behind the scenes’. This makes the tool very user friendly and suitable for todays tutorial. There is extensive documentation on the use of the pipeline. The following tutorial is the closest to what we will do today: https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

Set working directory

You may need to run this on the terminal, since in RStudio, the working directory is set to your home directory after the chunk is run! An alternative is to use the pull down menu at the top of the bottom-right window on RStudio. Click Files/More to see the option

setwd("/Users/obasak/Documents/1_projects/Teaching/BMW_Bioinformatics/2026/COO/")

Load libraries

library(dplyr,verbose = FALSE)
library(ggplot2,verbose = FALSE)
library(stringr,verbose = FALSE)
library(Seurat,verbose = FALSE)
library(RColorBrewer,verbose = FALSE)
library(useful,verbose = FALSE)
library(readr,verbose = FALSE)
library(hdf5r,verbose = FALSE)

2. Seurat analysis

For today’s tutorial, we will use the scRNAseq atlas of the adolescent mouse cortex published by Saunders et al 2018.

Aim

Process the scRNAseq data of the adolescent mouse cortex in order to reveal major cell clusters and identify cell types

\(~\)
Why? A major use of the scRNAseq is cell type identification. To achieve this, you need to perform quality control steps and cluster analysis. To visualize the data, you will perform dimensionality reduction. Finally, you can plot marker genes that you find from the literature to reveal the cell type identity of clusters \(~\)

We will use the Seurat object that we have uploaded. This object made specially for the Seurat pipeline has a lot of ‘slots’ where information can be stored, and the architecture is a bit complicated. You do not need it for this tutorial, except what is mentioned

2.1 The dataset

The data (‘Linnerson_cortex_10X22.rds’) downloaded and processed into a ‘Seurat object’ to prevent technical errors caused by the ‘loom’ file format in several computers.

You can download them from Brightspace.

2.1.1 Load the dataset

Load the Seurat object that was saved as rds.

dataset <- readRDS(file = 'data/Linnerson_cortex_10X22.rds')
dataset

Note: The object contains data from 6658 cells (samples) and 27998 features (genes). There is 1 assay (RNA).

2.1.2 Check the metadata

This is where ‘cell level’ information is stored. This means there is one value for each cell.
We will look at the first 100 to get an idea what is in.

head(dataset@meta.data$Class, 100) 
head(dataset@meta.data[,1:6])

2.2 Quality metrics

** We skip this part for time reasons **

2.3 Filter low quality cells

** We skip this part for time reasons **

2.4 Normalise

In Seurat, standard preprocessing workflow is replaced by a single command. However, it is good to see this part to learn each step.

First, we will normalize the data. This is to get rid of the differences in total RNA counts between cells. In other words, we will equalize the total count number in each cell to a fixed number (e.g. 10000 RNA molecules per cell).

dataset <- NormalizeData(object = dataset, normalization.method = "LogNormalize", scale.factor = 10000)

2.5 Detection of variable genes across the single cells

We want to find out ‘informative genes’ that will explain biological differences to use in some of hte downstream applications. If a gene is expressed everywhere, that doesnt tell us much. However, if a gene is expressed in a subset of cells, this will cause ‘variation’.

We can detect these genes using FindVariableFeatures() function.

## Here, we select the top 1,000 highly variable genes (Hvg) for downstream analysis.
dataset <- FindVariableFeatures(object = dataset, selection.method = 'mean.var.plot', mean.cutoff = c(0.0125, 3), dispersion.cutoff = c(0.5, Inf))
length(x = VariableFeatures(object = dataset)) #3084

Discussion point 1

Go through the code: What is happening?
\(~\)

Identify the 10 most highly variable genes.

top10 <- head(VariableFeatures(dataset), 10)
top10

Now visualise.

## Plot
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(dataset)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot2

Note:Dispersion indicates variation, while red color shows ‘significantly variable’ genes

Select the top 1000 highly variable genes. For downstream applications, we will use these.

hv.genes <- head(rownames(HVFInfo(object = dataset)), 1000)

Discussion point 2

Check the dispersion table: What do you see?

What do we mean by variation?
\(~\)

2.6 Scale the data and get rid of the confounders -

We will use the ScaleData() function of hte Seurat package. The confounding factors can be discarded using ‘vars.to.regress’.

dataset <- ScaleData(object = dataset, features = hv.genes, vars.to.regress = c("percent.mt","nFeature_RNA"))

Dimentionality reduction

2.7 PCA analysis

2.7.1 Perform PCA

Performing Dimensionality reduction in Seurat is very simple!

dataset <- RunPCA(object = dataset, features = hv.genes, verbose = FALSE,npcs = 50)

2.7.2 Plot PCA results

plot1 <- DimPlot(object = dataset, reduction = 'pca',dims = c(1,2))
plot1

It is also easy to find out genes that contribute to the + and - direction of the top PCs! We can also plot the level of contribution.

print(dataset[["pca"]], dims = 1:5, nfeatures = 5) # First box
VizDimLoadings(dataset, dims = 1:2, reduction = "pca") # Second box

Discussion point 3

What can we learn from the PC1 versus PC2 plot?
What do the genes represent?
What can we learn from by looking at these genes?
\(~\)

2.7.3 Visualise how genes contribute to PCs on heatmaps

We can use an integrated function of Seurat to plot heatmaps to visualize genes that drive different principle components.

Plot some of these genes on PCA plots to see how their expression is distributed along different PCs. For this, we will look at the first 4 PCs, indentify genes that are highest on + and - axis and plot them on PC1/PC2 and PC3/PC4 plots.

The first two components:

PCHeatmap(dataset,dims = c(1,2),ncol =2)

Cldn5 and Mog genes look interesting:

plot_f1 <- FeaturePlot(dataset,features = c('Cldn5','Mog'),reduction = 'pca',dims = c(1,2),cols = c('gray','blue','red'))
plot_f2 <- FeaturePlot(dataset,features = c('Cldn5','Mog'),reduction = 'pca',dims = c(3,4),cols = c('gray','blue','red'))
plot_f1 / plot_f2

Now components 5 and 8:

PCHeatmap(dataset,dims = c(5,8),ncol =2)

Pprm and Tmem119 genes look interesting:

plot_f3 <- FeaturePlot(dataset,features = c('Rprm','Tmem119'),reduction = 'pca',dims = c(5,8),cols = c('gray','blue','red'))
plot_f4 <- FeaturePlot(dataset,features = c('Rprm','Tmem119'),reduction = 'pca',dims = c(1,2),cols = c('gray','blue','red'))
plot_f3 / plot_f4

2.8 Find clusters using SNN Graph Construction and the louvain algorithm

We will use a graph-based clustering algorithm discussed at the lecture.

2.8.1 SNN Graph Construction

We need to build a neighborhood graph. In this network graph, each cell will be a node, and their similarity in high dimensional space will become their edges.

One could say that cells closest to each other reside in a neighborhood.

dataset <- FindNeighbors(object = dataset, dims = 1:20)

2.8.2 Find clusters using the louvain algorithm

dataset <- FindClusters(object = dataset, resolution = 0.6) # changing the resolution will change the size/number of clusters! we can try the following; c(0.6, 0.8, 1)
# we can plot the number of genes per cluseter using the violing plot
VlnPlot(object = dataset, features = "nCount_RNA")

2.9 UMAP and t-SNE analysis

2.9.1 Calculate the coordinates

We will use the top PCs to calculate the umap and tSNE coordinates. You can change the numnber of PCs based on your observation in Practice 2.

dataset <- RunUMAP(object = dataset, reduction = "pca", dims = 1:20, verbose = FALSE)
dataset <- RunTSNE(object = dataset, dims = 1:20, verbose = FALSE)

2.9.2 Plot the results

Adjust the DimPlot function to show the PCA, UMAP and t-SNE results

DimPlot(object = dataset, reduction = 'umap')
DimPlot(object = dataset, reduction = 'pca')
DimPlot(object = dataset, reduction = 'tsne')

Discussion point 4

Looking at the clusters that you have identified, what differences do you see between the dimentionality reduction methods? Do you have a preference to use one instead of the other? Why?
\(~\)

3. PRACTICE!

Practice

Identify cell types!

Check Use Google/ChatGPT to find marker genes for neurons, inhibitory neurons, astrocytes and oligodendrocytes.

Hint: The following papers have marker genes for relevant cell types:

Plot the expression of these marker genes using the VlnPLot() function.

** What do you see? Can you identify different cell types clearly, or is there anything missing? ** \(~\)

END OF THE PRACTICAL!

The following parts are extras. If you have time, please follow them :::

4. Example questions

The following is a question from teh last year. We didn’t go into details of k-Means this year, thus the question may look out of context. But it will give you an idea of what can be expected. A indicated the answer

Question 1

Please list four important facts about the k-means algorithm on the following topics: i) What is it used for? > A: For classification of data into groups/clusters

  1. Please explain the important steps > A: Determine the expected number of clusters. Take random points in the data and calculate the distance of each point to these random points to assign clusters. Then calculate the center of the cluster. Finally, repeat this process until there is no change in the centeral point, meaning that a stability is reached.

  2. Name one of the drawbacks > A: needs estimation of the number of clusters beforehand. Cannot work on all different ‘shapes’ of data. Cannot find outliers. Does not show how similar each cluster is.

  3. If you run the k-means algorithm on the same data two different times, do you always get the same results? Why/why not? > A: No. The process is stochastic, starts randomly and is repeated many times until stability is reached. The final results will, in most cases, be different.


Here is another question for you to think about:

Question 2

Which steps of the single cell RNA sequencing analysis aims at getting rid of the differences in number of reads between cells (e.g. a cell with 1000 features and another with 5000 features)?

  1. Scaling
  2. Normalization
  3. Regression
  4. Dimensionality reduction
  5. Clustering

A: I wont provide the answer for this one