Introductory Bioinformatics for Neuroscience
MSI Research Informatics
2026-08-04
1 Set up
To analyze single nuclei/cell data you will most likely need
specialized software (aka. “packages”). We have already installed what
you need for this tutorial, so you just need to tell R where to find it
and “load” it. The .libPaths function helps you see or
change the locations where R will look for software, we will use it to
include a new location stored in the tutorial_packages
variable, then we will load the packages we will use:
tutorial_packages <- "/common/tutorials/bioinformatics/single_cell/software/r-rich4.6_lib/"
.libPaths(c(tutorial_packages, .libPaths()))
.libPaths()## [1] "/common/tutorials/bioinformatics/single_cell/software/r-rich4.6_lib"
## [2] "/users/0/calix015/R/x86_64-pc-linux-gnu-library/4.6"
## [3] "/opt/R-4.6.0/lib64/R/library"library(ggplot2)
library(tidyverse)
library(Seurat)When you generate results, they will be stored by default in your
current working directory, so it is handy to learn how to find or change
it. This uses the getwd and setwd
functions.
# Here is where your outputs will be stored by default:
getwd()## [1] "/users/0/calix015"# If you want to change it:
project_dir <- "/projects/standard/riss/calix015/sc_tutorial/"
setwd(project_dir)Finally, for this tutorial we will be using data downloaded from GEO
with the accession GSE247974. Note that we are assigning the location of
this data (aka. its “path”) to the variable data_dir so
that we can reuse it easily.
data_dir <- "/common/tutorials/bioinformatics/single_cell/data/neuro_dataset/"
# You can see the contents of this directory with the "list.files" function:
list.files(data_dir)## [1] "GSE247974_D3" "GSE247974_D3.Rds"list.files(paste0(data_dir,"GSE247974_D3"))## [1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz"2 How to load the count data?
This depends on the format of the input. The Seurat
package has different functions for different types of data and you can
browse and find more information here. For this data
set (10X) we will use the Read10X function to store the
data in the so variable:
count_mtx <- Read10X(data.dir = paste0(data_dir,"GSE247974_D3"))We can see details about how to use this function (or any other) and
the default parameters with the ? command:
?Read10XNow we will inspect the size and contents of the count matrix we just
loaded. We will use a subset of the rows and columns to get a little
preview instead of showing the entire matrix with
[ , ]:
dim(count_mtx)## [1] 36601 4245count_mtx[1:10,1:5]## 10 x 5 sparse Matrix of class "dgCMatrix"
## AAACCCAAGAGGGTAA-1 AAACCCACACAACCGC-1 AAACCCAGTAAGCTCT-1
## MIR1302-2HG . . .
## FAM138A . . .
## OR4F5 . . .
## AL627309.1 . . .
## AL627309.3 . . .
## AL627309.2 . . .
## AL627309.5 . . .
## AL627309.4 . . .
## AP006222.2 . . .
## AL732372.1 . . .
## AAACCCATCCGTGGGT-1 AAACCCATCGCCATAA-1
## MIR1302-2HG . .
## FAM138A . .
## OR4F5 . .
## AL627309.1 . .
## AL627309.3 . .
## AL627309.2 . .
## AL627309.5 . .
## AL627309.4 . .
## AP006222.2 . .
## AL732372.1 . .This small example we have chosen for the tutorial has information
from 36,601 genes and 4,245 barcodes. Note that most of the values are
., which represent zeros. This is what a “sparse” data set
looks like. Also note that the row names are gene names and the column
names are barcodes.
3 Preprocessing overview
3.1 Format the counts matrix
To analyze a data set like this we should first format it so that the
packages we will use can access the information. The
CreateSeuratObject function does that:
so <- CreateSeuratObject(count_mtx)
so## An object of class Seurat
## 36601 features across 4245 samples within 1 assay
## Active assay: RNA (36601 features, 0 variable features)
## 1 layer present: countsSince we know there are a lot of zeros, we can use some minimum
thresholds to filter out the barcodes that are most likely empty
droplets. We will rewrite the so variable:
# Create a SO with filters
so <- CreateSeuratObject(counts = count_mtx, min.features = 10, min.cells = 5)
so## An object of class Seurat
## 20759 features across 4245 samples within 1 assay
## Active assay: RNA (20759 features, 0 variable features)
## 1 layer present: counts3.2 Structure of a Seurat object (so)
You can browse the structure of the data by clicking on its name in
the Data panel or with the View function. The
head function shows you a little preview instead of the
entirety of the data. We will peek into a few key elements:
# all the genes
rownames(so) %>% head()## [1] "AL627309.5" "LINC01409" "LINC01128" "LINC00115" "FAM41C"
## [6] "AL645608.6"# the counts
so@assays$RNA$counts[1:10,1:5]## 10 x 5 sparse Matrix of class "dgCMatrix"
## AAACCCAAGAGGGTAA-1 AAACCCACACAACCGC-1 AAACCCAGTAAGCTCT-1
## AL627309.5 . . .
## LINC01409 . . .
## LINC01128 . . 1
## LINC00115 . . .
## FAM41C . . .
## AL645608.6 . . .
## AL645608.2 . . .
## AL645608.4 . . .
## LINC02593 . . .
## SAMD11 1 . .
## AAACCCATCCGTGGGT-1 AAACCCATCGCCATAA-1
## AL627309.5 . .
## LINC01409 . .
## LINC01128 . .
## LINC00115 . .
## FAM41C . .
## AL645608.6 . .
## AL645608.2 . .
## AL645608.4 . .
## LINC02593 . .
## SAMD11 . .# the metadata
so@meta.data %>% head()3.3 QC to clustering
Please check out the Seurat tutorials for detailed information. Briefly, this process involves a quality control check and filtering out barcodes that don’t contain acceptable cells/nuclei, count normalization, exploring the variability of the data, scaling, dimensionality reduction, clustering, cluster identification and sometimes sample integration.
A template workflow would look like this:
# calculate filtering metrics:
so[["percent.mt"]] <- PercentageFeatureSet(so, pattern = "^MT-")
# plot values and select thresholds:
VlnPlot(so,c("nCount_RNA", "nFeature_RNA","percent.mt"))
# filter:
so <- subset(so,
subset = nFeature_RNA > 1000 & nFeature_RNA < 7500 & nCount_RNA < 50000)
# Normalize:
so <- NormalizeData(so)
# Inspect variability:
so <- FindVariableFeatures(so)
# Scale:
so <- ScaleData(so, features = rownames(so))
# Reduce dimensionality:
so <- RunPCA(so,
features = VariableFeatures(object = so))
# Define dimensionality threshold:
ElbowPlot(so)
ndims <- 15
# Cluster:
so <- FindNeighbors(so, dims = 1:ndims)
so <- FindClusters(so, resolution = 0.5)
# Calculate other data projections for plotting:
so <- RunUMAP(so, dims = 1:ndims)But be aware that the thresholds and options chosen in these functions should mold to your data’s characteristics. Please think carefully about these steps because they are the base for the results you will get.
Finally, for cluster or cell identification we need information about the samples and the expected cell types, marker genes or a really good database to compare the data against. This process is quite particular to your samples, iterative and sometimes quite difficult. We talk a little more in depth about it in our “Introduction to Single Cell Genomics” tutorial.
4 Pre-analyzed object
In this part we will show you how to open a pre-analyzed object in case you need to, for example, generate plots of your favorite genes or cell types.
# load the pre-processed and labeled data
so_done <- readRDS(paste0(data_dir,"GSE247974_D3.Rds"))
so_done## An object of class Seurat
## 43120 features across 3908 samples within 2 assays
## Active assay: SCT (20675 features, 3000 variable features)
## 3 layers present: counts, data, scale.data
## 1 other assay present: RNA
## 2 dimensional reductions calculated: pca, umapAfter analysis we have more layers with gene expression information, a set of variable genes, dimensional reductions and more information in the metadata:
so_done@meta.data %>%head()The object meta data is nothing more than a data frame with the cell barcodes as row names and different cell attributes in the columns, so we can count and summarize this information as needed. For example, we can look at the “celltype” column and count how many cells of each type are present:
table(so_done$celltype) ##
## BG Brain-Stem Choroid
## 126 46 522
## Committed-OPC Div-Choroid Div-VZ
## 4 256 73
## eCN/Unibrush Endo GC
## 49 481 86
## GCP Glia iCN
## 190 161 3
## immature-iCN/immature-PC Meninges Micro
## 137 104 8
## MLI PC PIP
## 1 35 35
## Progenitors RL Roof-Plate
## 748 47 310
## VZ
## 32To assign each cell an identity value stored in a metadata column, we
use the Idents function:
Idents(so_done) <- so_done$celltype4.1 Subset cells
You can also use the metadata to subset a group of cells that most interest you. For example, let’s subset all cell types that have at least 10 representatives:
# select the ids that you want with a condition:
ids_ok <- table(so_done$celltype) > 10
ids_ok <- names(ids_ok)[ids_ok]
#subset:
so_done_subset <- subset(so_done, idents = ids_ok)4.2 Plot cell projections
The fisrt 2 dimensions after the dimensionality reduction step are usually chosen for a 2D representation of your data. Commonly, a PCA reduction is the basis for all downstream analysis, let’s see how the cells look when projected in a 2D space:
DimPlot(so_done_subset, reduction = "pca")
Often, you’ll see another very common projection strategy called UMAP:
DimPlot(so_done_subset)
Note that neither of these plots tells you the full story of the data,
it’s just a “summary picture” and a very common mistake is to
misinterpret the distance between clusters in a UMAP with the difference
between cells. This is generally not correct and can lead you to wrong
conclusions.
You can customize your plot changing the function parameters (use
?DimPlot). Here I am removing the legend and annotating the
cell type over the plot:
DimPlot(so_done_subset,label = TRUE, pt.size = 0.5) + NoLegend()
4.3 Plot gene expression
Two very common plotting options to summarize the gene expression of your genes of interest across different samples or cell groups are the violin and dot plots. Let’s make some with the most variable genes in this data set:
# top 5 var genes
top5 <- so_done_subset@assays$SCT@var.features[1:5]
# Violin of the most variable:
VlnPlot(so_done_subset, top5[1])+NoLegend()+xlab("")
# Dotplot:
DotPlot(so_done_subset, top5)
Let’s improve this dotplot a little:
holder <- DotPlot(so_done_subset, top5)
holder <- holder+ggtitle("customizing my plot")+
scale_color_gradient(low = "grey",high = "red")+
xlab("")+
ylab("")+
theme(axis.text = element_text(size=8))
holder
If you’d like to see a particular gene’s expression data or a metadata column with quantitative information projected on the UMAP plot, you can use a “feature” plot:
colnames(so_done_subset@meta.data)## [1] "orig.ident" "nCount_RNA" "nFeature_RNA" "percent.mt"
## [5] "nCount_SCT" "nFeature_SCT" "SCT_snn_res.0.5" "seurat_clusters"
## [9] "celltype"# Gene expression:
FeaturePlot(so_done_subset, top5[1])
# Metadata column:
FeaturePlot(so_done_subset, "nCount_SCT")
5 Save your work!
Remember to save your results and plots before your session ends.
# Save the SO:
saveRDS(so_done_subset, paste0(project_dir,"subset.rds"))
# Save your custom plots:
ggsave(paste0(project_dir,"custom_plot.png"), plot = holder)