Clustering 3K PBMCs with Bioconductor

Author

Kevin Rue-Albrecht

Generated: 2026-07-17 11:15:29

Context

This tutorial demonstrates a single-cell analysis workflow using Bioconductor packages, while discussing its implementation as a Galaxy Training Network (GTN) tutorial.

Libraries

library(bluster)
library(cowplot)
library(DropletUtils)
library(LoomExperiment)
library(scater)
library(scran)
library(scuttle)

Data

Data download

To create a parallel tutorial to the existing tutorial Clustering 3K PBMCs with Seurat, input data files were downloaded from these Zenodo links:

https://zenodo.org/record/3581213/files/genes.tsv
https://zenodo.org/record/3581213/files/barcodes.tsv
https://zenodo.org/record/3581213/files/matrix.mtx

if (!dir.exists("data")) {
    dir.create("data")
}
if (!file.exists("data/genes.tsv")) {
    download.file("https://zenodo.org/record/3581213/files/genes.tsv", "data/genes.tsv")
}
if (!file.exists("data/barcodes.tsv")) {
    download.file("https://zenodo.org/record/3581213/files/barcodes.tsv", "data/barcodes.tsv")
}
if (!file.exists("data/matrix.mtx")) {
    download.file("https://zenodo.org/record/3581213/files/matrix.mtx", "data/matrix.mtx")
}

Create a SingleCellExperiment Object

The common data structure adopted by Bioconductor packages to store the matrix as well as gene and cell annotations is called SingleCellExperiment.

We can import the matrix and annotations of genes and cells into an SingleCellExperiment object using the function read10xCounts() of the package DropletUtils.

We set the samples= argument to the directory that contains the three 10x files. Changing row.names= from the default "id" to "symbol" ensures that human-readable gene symbols are used instead of the database-friendlier Ensembl gene identifiers. Either way, both Ensembl gene identifiers and gene symbols will be stored in the rowData() component of the SingleCellExperiment and available at all times.

In the context of a Galaxy workflow, we want to save the R object to a file that can be re-used by downstream tasks. R objects are traditionally saved to .RData or .rds files. However those file formats are specific to R and deemed insecure.

To be saved to the safer loom file format, SingleCellExperiment objects must be converted to the SingleCellLoomExperiment class before they can be saved to disk using the export() function from the BiocIO package.

Altogether, the first tool in our Galaxy workflow should execute the following code:

library(DropletUtils)
library(LoomExperiment)
sce <- DropletUtils::read10xCounts(
    samples = "data/",
    row.names = "symbol"
)
if (!dir.exists("outputs")) {
dir.create("outputs")
}
scle <- as(sce, "SingleCellLoomExperiment")
if (!file.exists("outputs/sce.loom")) {
    export(object = scle, con = "outputs/sce.loom", format = "loom")
}
rm(list = ls())

For this task, the tool DropletUtils Read10x can be used.

However:

This tool produces an rdata file. RData objects are deemed insecure as discussed in this GitHub issue. The tool should be updated to produce a loom file containing a LoomExperiment object.
This tool does not offer the option to use gene symbols as rownames. The default Ensembl gene identifiers are not immediately interpretable and complicate the interpretation of results later in the workflow (e.g., list of highly variable genes, cluster marker genes). The tool should be updated to offer the possibility of using gene symbols as rownames for the SingleCellExperiment. It is worth pointing out that both Ensembl gene identifiers and gene symbols are stored in the rowData component of the SingleCellExperiment object, meaning they remain available throughout the analysis

Inspect the `SingleCellExperiment` object.

Given that R sessions do not persist between steps of a Galaxy workflow, a Galaxy tool would need to execute the following code to re-import the SingleCellLoomExperiment from the file written in the previous step before displaying a summary view of it.

In the context of a Galaxy workflow, the output would then be saved to a .txt file that the user could inspect after the job completes.

library(LoomExperiment)
sce <- import('outputs/sce.loom', format = "loom", type = "SingleCellLoomExperiment")
capture.output(print(sce), file = "outputs/sce_summary.txt")
rm(list = ls())

For this task, a tool similar to the Seurat Data Management tool method “Inspect Seurat Object” is needed.

Preprocessing

Computation of QC metrics

Identify mitochondrial genes.

In preparation for the calculation of quality control metrics, we need to identify mitochondrial genes in our dataset.

Conveniently, human genes that are encoded in the mitochondrial DNA (rather than in the cell nucleus) have names beginning with MT-.

The following code identifies all the rownames that match this pattern in our dataset and writes them to a .txt file.

library(LoomExperiment)
sce <- import('outputs/sce.loom', format = "loom", type = "SingleCellLoomExperiment")
mt_gene_ids <- rownames(sce)[grep('^MT-', rownames(sce))]
cat(mt_gene_ids, file = "outputs/mt_gene_ids.txt", sep = "\n")
cat(mt_gene_ids, sep = "\n")

MT-ND1
MT-ND2
MT-CO1
MT-CO2
MT-ATP8
MT-ATP6
MT-CO3
MT-ND3
MT-ND4L
MT-ND4
MT-ND5
MT-ND6
MT-CYB

rm(list = ls())

Calculate the Proportion of Mitochondrial Reads

Equipped with a SingleCellLoomExperiment object and a file that contains the list of rownames that correspond to the mitochondrial genes, we can use the addPerCellQCMetrics() function from the scuttle package to compute some generic quality control metrics alongside metrics focusing on mitochondrial genes.

In particular, the percentage of UMI assigned to mitochondrial genes is frequently used as a measure of cell integrity allowing us to remove droplets that contain damaged or otherwise suspicious cells.

library(LoomExperiment)
library(scuttle)
sce <- import('outputs/sce.loom', format = "loom", type = "SingleCellLoomExperiment")
mt_gene_ids <- scan("outputs/mt_gene_ids.txt", what = "character", quiet = TRUE)
gene_subsets <- list()
gene_subsets[["MT"]] <- mt_gene_ids
sce <- scuttle::addPerCellQCMetrics(
    x = sce,
    subsets = gene_subsets
)
if (!file.exists("outputs/sce_after_qc.loom")) {
    export(object = sce, con = "outputs/sce_after_qc.loom", format = "loom")
}
rm(list = ls())

For this task, the tool Scater: calculate QC metrics exists, but:

It requires the expression matrix in tabular format. The tool should be updated to take a loom file containing a SingleCellLoomExperiment object.

Notes:

Single-cell quality control with scater tutorial provides a precomputed file of mitochondrial genes.

Inspect the QC metrics

The newly computed per-cell quality control metrics are stored in the colData component of the object. This a tabular DataFrame structure which can be inspected by re-importing the SingleCellLoomExperiment object and printing the colData() component.

In the context of a Galaxy workflow, the output would then be saved to a .txt file that the user could inspect after the job completes.

library(LoomExperiment)
sce <- import('outputs/sce_after_qc.loom', format = "loom", type = "SingleCellLoomExperiment")
write.table(
    as.data.frame(colData(sce)), file = "outputs/sce_coldata.tsv", sep = "\t", 
    quote = FALSE, row.names = FALSE
)
rm(list = ls())

Filtering of low-quality cells

Visualise QC Metrics.

Having computed the quality control metrics and stored them in the colData component of the SingleCellLoomExperiment object, we can reload that object and use the plotColData() function of the scater package to visualise the quality control metrics to determine appropriate thresholds for filtering out low-quality cells.

Which metrics are available can be determined by inspecting the colData() component of the object as shown above.

Given that plotColData() can only plot one QC metric at a time on the y-axis, we use the cowplot package to combine them into a single figure.

library(LoomExperiment)
library(scater)
sce <- import('outputs/sce_after_qc.loom', format = "loom", type = "SingleCellLoomExperiment")
plot_list <- list()
plot_list[["detected"]] <- scater::plotColData(
    object = sce,
    y = "detected"
)
plot_list[["sum"]] <- scater::plotColData(
    object = sce,
    y = "sum"
)
plot_list[["subsets_MT_percent"]] <- scater::plotColData(
    object = sce,
    y = "subsets_MT_percent"
)
cowplot::plot_grid(plotlist = plot_list, nrow = 1)

rm(list = ls())

The plotColData() function can also plot one QC metric on the y-axis against another QC metric on the x-axis. This can be additional insights into the relationship between QC metrics and help determine appropriate thresholds for filtering out low-quality cells.

library(LoomExperiment)
library(scater)
sce <- import('outputs/sce_after_qc.loom', format = "loom", type = "SingleCellLoomExperiment")
plot_list <- list()
plot_list[["detected"]] <- scater::plotColData(
    object = sce,
    y = "detected",
    x = "sum"
)
plot_list[["subsets_MT_percent"]] <- scater::plotColData(
    object = sce,
    y = "subsets_MT_percent",
    x = "sum"
)
cowplot::plot_grid(plotlist = plot_list, nrow = 1)

rm(list = ls())

Notes:

In contrast to the ‘violin plot’ use case above, where users can give a list of column names that the wrapper can simply loop over, this ‘scatter plot’ use case requires two features per plot - one of the y-axis and one for the x-axis - which immediately complicates the idea of producing multiple plots for multiple combinations of features on the X-Y axes.
Again, taking inspiration from the Seurat Visualize, the tool should only allow users to specify one feature for the x-axis and one feature for the y-axis, forcing them to call scater::plotColData() once per combination of features on the X-Y axes.
The code demonstrated above uses cowplot::plot_grid() to combine multiple plots into a single figure only for the purpose of producing an image to paste in the GTN tutorial. In a Galaxy workflow, the tool should produce one image per call to scater::plotColData() and not attempt to combine them into a single figure. {: .comment}

Filter Out Low Quality Cells

Having visualised the quality control metrics, we can now filter out low-quality cells based on thresholds for the number of detected genes and the percentage of mitochondrial reads.

library(LoomExperiment)
sce <- import('outputs/sce_after_qc.loom', format = "loom", type = "SingleCellLoomExperiment")
keep_cells <- sce$detected >= 200 & sce$detected <= 2500 & sce$subsets_MT_percent <= 5
sce <- sce[, keep_cells]
if (!file.exists("outputs/sce_after_qc_filter.loom")) {
    export(object = sce, con = "outputs/sce_after_qc_filter.loom", format = "loom")
}
rm(sce)

Notes:

Not every user will want to apply the same types of filters. Typically there are three metrics (total UMI, number of genes detected, mitochondrial content). For each metrics, users might choose to apply a lower cutoff and a higher cutoff. The Galaxy tool should offer the flexibility to filter on any of those six metrics. Take inspiration from the equivalent Seurat tool that already exists.

Visualise QC Metrics After Filtering

We can then produce the same plots as earlier to visualise the new distribution of quality control metrics in our data.

library(LoomExperiment)
library(scater)
sce <- import('outputs/sce_after_qc_filter.loom', format = "loom", type = "SingleCellLoomExperiment")
plot_list <- list()
plot_list[["detected"]] <- scater::plotColData(
    object = sce,
    y = "detected"
)
plot_list[["sum"]] <- scater::plotColData(
    object = sce,
    y = "sum"
)
plot_list[["subsets_MT_percent"]] <- scater::plotColData(
    object = sce,
    y = "subsets_MT_percent"
)
cowplot::plot_grid(plotlist = plot_list, nrow = 1)

rm(list = ls())

library(LoomExperiment)
library(scater)
sce <- import('outputs/sce_after_qc_filter.loom', format = "loom", type = "SingleCellLoomExperiment")
plot_list <- list()
plot_list[["detected"]] <- scater::plotColData(
    object = sce,
    y = "detected",
    x = "sum"
)
plot_list[["subsets_MT_percent"]] <- scater::plotColData(
    object = sce,
    y = "subsets_MT_percent",
    x = "sum"
)
cowplot::plot_grid(plotlist = plot_list, nrow = 1)

rm(list = ls())

Further Preprocessing

Log-normalisation

Having filtered droplets that do not pass our chosen quality control filters, the next step in a standard single-cell RNA-sequencing analysis workflow is normalisation.

In Bioconductor, the simple log-normalisation method is implemented in the function logNormCounts() of the package scuttle.

library(LoomExperiment)
library(scuttle)
sce <- import('outputs/sce_after_qc_filter.loom', format = "loom", type = "SingleCellLoomExperiment")
sce <- scuttle::logNormCounts(x = sce)
dimnames(sce@assays@data$counts) <- list(NULL, NULL)
dimnames(sce@assays@data$logcounts) <- list(NULL, NULL)
if (!file.exists("outputs/sce_after_lognorm.loom")) {
    export(object = sce, con = "outputs/sce_after_lognorm.loom", format = "loom")
}

Model gene variance

We can then model gene variance using the matrix normalised of normalised counts using the modelGeneVar() function of the scran package.

library(LoomExperiment)
library(scuttle)
sce <- import('outputs/sce_after_lognorm.loom', format = "loom", type = "SingleCellLoomExperiment")
allf <- scran::modelGeneVar(x = sce, assay.type = "logcounts")
allf_df <- as.data.frame(allf)
allf_df <- data.frame(gene = rownames(allf_df), allf_df, row.names = NULL, check.names = FALSE)
write.table(
    allf_df, file = "outputs/scuttle_model_gene_var.tsv", sep = "\t",
    quote = FALSE, row.names = FALSE, col.names = TRUE
)

hvgs <- scran::getTopHVGs(allf, n = 2000)

message("No direct equivalent to Seurat ScaleData")

No direct equivalent to Seurat ScaleData

plot(allf$mean, allf$total)
curve(metadata(allf)$trend(x), add=TRUE, col="dodgerblue")

Dimensionality Reduction

Principal Component Analysis

Perform the PCA

sce <- scater::runPCA(
    x = sce,
    ncomponents = 50,
    subset_row = hvgs
)
sce

class: SingleCellLoomExperiment 
dim: 32738 2638 
metadata(5): CreatedWith LOOM_SPEC_VERSION LoomExperiment-class
  MatrixName Samples
assays(2): counts logcounts
rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
rowData names(3): ID Symbol subsets_MT
colnames: NULL
colData names(9): Barcode Sample ... sum total
reducedDimNames(1): PCA
mainExpName: NULL
altExpNames(0):
rowGraphs(0): NULL
colGraphs(0): NULL

Visualise the PCA Results

Visualise the PCA Results - Dimensional Loadings

message("No direct equivalent to Seurat::VizDimLoadings")

No direct equivalent to Seurat::VizDimLoadings

Visualise the PCA Results - DimPlot

scater::plotPCA(sce)

Visualise the PCA Results - Feature Plots

plot_list <- list()
for (gene_id in c("CST3", "CD79A", "HLA-DQA1", "MALAT1", "NKG7", "PPBP")) {
    plot_list[[gene_id]] <- scater::plotPCA(
        object = sce,
        ncomponents = c(1, 2),
        colour_by = gene_id
    )
}
cowplot::plot_grid(plotlist = plot_list, nrow = 2)

Same for PC2 and PC3.

plot_list <- list()
for (gene_id in c("CST3", "CD79A", "HLA-DQA1", "MALAT1", "NKG7", "PPBP")) {
    plot_list[[gene_id]] <- scater::plotPCA(
        object = sce,
        ncomponents = c(2, 3),
        colour_by = gene_id
    )
}
cowplot::plot_grid(plotlist = plot_list, nrow = 2)

Visualise the PCA Results - Heatmaps

message("No direct equivalent to Seurat::DimHeatmap")

No direct equivalent to Seurat::DimHeatmap

Determine the ‘dimensionality’ of the dataset

Make an Elbow Plot

percent.var <- attr(reducedDim(sce), "percentVar")
plot(percent.var, log="y", xlab="PC", ylab="Variance explained (%)")

Cluster the Cells

colData(sce)[["cluster"]] <- scran::clusterCells(
    x = sce,
    use.dimred = "PCA",
    BLUSPARAM = bluster::SNNGraphParam(
        k = 10,
        type="rank",
        cluster.fun = "louvain",
        cluster.args = list(
            resolution = 0.7
        )
    )
)
sce

class: SingleCellLoomExperiment 
dim: 32738 2638 
metadata(5): CreatedWith LOOM_SPEC_VERSION LoomExperiment-class
  MatrixName Samples
assays(2): counts logcounts
rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
rowData names(3): ID Symbol subsets_MT
colnames: NULL
colData names(10): Barcode Sample ... total cluster
reducedDimNames(1): PCA
mainExpName: NULL
altExpNames(0):
rowGraphs(0): NULL
colGraphs(0): NULL

Run non-linear dimensional reduction (UMAP/tSNE)

Perform Dimensional Reduction with UMAP

sce <- scater::runUMAP(
    x = sce,
    dimred = "PCA",
    pca = 20
)
sce

class: SingleCellLoomExperiment 
dim: 32738 2638 
metadata(5): CreatedWith LOOM_SPEC_VERSION LoomExperiment-class
  MatrixName Samples
assays(2): counts logcounts
rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
rowData names(3): ID Symbol subsets_MT
colnames: NULL
colData names(10): Barcode Sample ... total cluster
reducedDimNames(2): PCA UMAP
mainExpName: NULL
altExpNames(0):
rowGraphs(0): NULL
colGraphs(0): NULL

Visualise UMAP

Visualise the UMAP

scater::plotUMAP(
    object = sce,
    colour_by = "cluster",
    text_by = "cluster"
)

plot_list <- list()
for (gene_id in c("CST3", "CD79A", "HLA-DQA1", "MALAT1", "NKG7", "PPBP")) {
    plot_list[[gene_id]] <- scater::plotUMAP(
        object = sce,
        colour_by = gene_id
    )
}
cowplot::plot_grid(plotlist = plot_list, nrow = 2)

Other ways to visualise clusters

Visualise Expression across Clusters with a Violin Plot

plotExpression(
    object = sce,
    features = c("CST3", "CD79A", "HLA-DQA1", "MALAT1", "NKG7", "PPBP"),
    ncol = 3,
    x = "cluster",
    colour_by = "cluster"
)

Find Marker Genes

Find Positive Markers for Every Cluster Compared to the Rest

Use `findMarkers` to Compare All Clusters

find_markers_results <- findMarkers(
    x = sce,
    groups = colData(sce)[["cluster"]],
    test.type = "wilcox",
    pval.type = "any"
)
find_markers_results

List of length 8
names(8): 1 2 3 4 5 6 7 8

could be converted to a base R list to safely write to RDS without triggering packages when re-importing in subsequent steps.

Find Markers of Cluster 2

Use `FindMarkers` on a Single Cluster

find_markers_results[["2"]]

DataFrame with 32738 rows and 11 columns
                 Top      p.value          FDR summary.AUC     AUC.1     AUC.3
           <integer>    <numeric>    <numeric>   <numeric> <numeric> <numeric>
FCGR3A             1  5.49380e-95  3.67053e-92   0.0178499  0.490248  0.447806
SDPR               1  2.81402e-78  1.09673e-75   0.0000000  0.495915  0.493204
GZMA               1  1.09260e-93  6.74900e-91   0.0332670  0.450028  0.484528
HLA-DRB1           1 7.01240e-162 1.14786e-157   0.9846097  0.984610  0.711375
CCL5               1 1.45422e-114 1.44267e-111   0.0108949  0.430930  0.490367
...              ...          ...          ...         ...       ...       ...
AC145205.1     32734            1            1         0.5       0.5       0.5
BAGE5          32735            1            1         0.5       0.5       0.5
CU459201.1     32736            1            1         0.5       0.5       0.5
AC002321.2     32737            1            1         0.5       0.5       0.5
AC002321.1     32738            1            1         0.5       0.5       0.5
               AUC.4     AUC.5     AUC.6     AUC.7     AUC.8
           <numeric> <numeric> <numeric> <numeric> <numeric>
FCGR3A     0.0973761 0.3974787 0.0178499  0.502627  0.518841
SDPR       0.4901961 0.4952381 0.4746835  0.495274  0.000000
GZMA       0.0332670 0.1129331 0.4745001  0.478266  0.505797
HLA-DRB1   0.9858293 0.9609110 0.7554119  0.988455  0.984980
CCL5       0.1828739 0.0108949 0.4962209  0.498019  0.000000
...              ...       ...       ...       ...       ...
AC145205.1       0.5       0.5       0.5       0.5       0.5
BAGE5            0.5       0.5       0.5       0.5       0.5
CU459201.1       0.5       0.5       0.5       0.5       0.5
AC002321.2       0.5       0.5       0.5       0.5       0.5
AC002321.1       0.5       0.5       0.5       0.5       0.5

[… skipped over a few variations for marker detection …]

Identify Cell Types

Make Violin Plots of B Cell Markers

plotExpression(
    object = sce,
    features = c("CD79A", "MS4A1"),
    ncol = 2,
    x = "cluster",
    colour_by = "cluster"
)

Make Violin Plots of T Cell Markers

plotExpression(
    object = sce,
    features = c("IL7R", "CCR7", "S100A4", "CD8A"),
    ncol = 4,
    x = "cluster",
    colour_by = "cluster"
)

Colour the UMAP Plot by Canonical Marker Expression

gene_ids <- c(
    "IL7R", "CCR7", "CD14", "LYZ", "S100A4", "MS4A1", "CD8A", "FCGR3A", "MS4A7",
    "GNLY", "NKG7", "FCER1A", "CST3", "PPBP"
)
plot_list <- list()
for (gene_id in gene_ids) {
    plot_list[[gene_id]] <- scater::plotUMAP(
        object = sce,
        colour_by = gene_id
    )
}
cowplot::plot_grid(plotlist = plot_list, nrow = 4)

Use Violin Plots to Compare Expression of Canonical Markers by Cluster

gene_ids <- c(
    "IL7R", "CCR7", "CD14", "LYZ", "S100A4", "MS4A1", "CD8A", "FCGR3A", "MS4A7",
    "GNLY", "NKG7", "FCER1A", "CST3", "PPBP"
)
plotExpression(
    object = sce,
    features = gene_ids,
    ncol = 4,
    x = "cluster",
    colour_by = "cluster"
)

Rename Clusters with Cell Types

new_cluster_labels <- c(
    "Memory CD4+ T",
    "B",
    "CD14+ Mono",
    "NK",
    "CD8+ T",
    "FCGR3A+ Mono",
    "Naive CD4+ T",
    "Platelet"
)
levels(colData(sce)[["cluster"]]) <- new_cluster_labels

Revisualise the UMAP with Cell Type Annotations

scater::plotUMAP(
    object = sce,
    colour_by = "cluster",
    text_by = "cluster"
)

Session information

sessionInfo()

R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.5.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/London
tzcode source: internal

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

other attached packages:
 [1] scran_1.38.1                scater_1.38.1              
 [3] ggplot2_4.0.3               scuttle_1.20.0             
 [5] LoomExperiment_1.28.0       BiocIO_1.20.0              
 [7] rhdf5_2.54.1                DropletUtils_1.30.0        
 [9] SingleCellExperiment_1.32.0 SummarizedExperiment_1.40.0
[11] Biobase_2.70.0              GenomicRanges_1.62.1       
[13] Seqinfo_1.0.0               IRanges_2.44.0             
[15] S4Vectors_0.48.1            BiocGenerics_0.56.0        
[17] generics_0.1.4              MatrixGenerics_1.22.0      
[19] matrixStats_1.5.0           cowplot_1.2.0              
[21] bluster_1.20.0             

loaded via a namespace (and not attached):
 [1] viridisLite_0.4.3         vipor_0.4.7              
 [3] farver_2.1.2              viridis_0.6.5            
 [5] R.utils_2.13.0            S7_0.2.2                 
 [7] fastmap_1.2.0             digest_0.6.39            
 [9] rsvd_1.0.5                lifecycle_1.0.5          
[11] cluster_2.1.8.2           statmod_1.5.2            
[13] magrittr_2.0.5            compiler_4.5.1           
[15] rlang_1.3.0               tools_4.5.1              
[17] igraph_2.3.3              yaml_2.3.12              
[19] FNN_1.1.4.1               knitr_1.51               
[21] labeling_0.4.3            S4Arrays_1.10.1          
[23] dqrng_0.4.1               DelayedArray_0.36.1      
[25] RColorBrewer_1.1-3        abind_1.4-8              
[27] BiocParallel_1.44.0       HDF5Array_1.38.0         
[29] withr_3.0.3               R.oo_1.27.1              
[31] grid_4.5.1                beachmat_2.26.0          
[33] Rhdf5lib_1.32.0           edgeR_4.8.2              
[35] scales_1.4.0              cli_3.6.6                
[37] rmarkdown_2.31            metapod_1.18.0           
[39] RSpectra_0.16-2           DelayedMatrixStats_1.32.0
[41] ggbeeswarm_0.7.3          stringr_1.6.0            
[43] parallel_4.5.1            BiocManager_1.30.27      
[45] XVector_0.50.0            vctrs_0.7.3              
[47] Matrix_1.7-5              jsonlite_2.0.0           
[49] BiocSingular_1.26.1       BiocNeighbors_2.4.0      
[51] ggrepel_0.9.8             irlba_2.3.7              
[53] beeswarm_0.4.0            h5mread_1.2.1            
[55] locfit_1.5-9.12           limma_3.66.0             
[57] glue_1.8.1                codetools_0.2-20         
[59] uwot_0.2.4                stringi_1.8.7            
[61] gtable_0.3.6              ScaledMatrix_1.18.0      
[63] htmltools_0.5.9           rhdf5filters_1.22.0      
[65] R6_2.6.1                  sparseMatrixStats_1.22.0 
[67] evaluate_1.0.5            lattice_0.22-9           
[69] R.methodsS3_1.8.2         renv_1.2.3               
[71] Rcpp_1.1.2                gridExtra_2.3.1          
[73] SparseArray_1.10.10       xfun_0.60                
[75] pkgconfig_2.0.3

Context

Libraries

Data

Data download

Create a SingleCellExperiment Object

Inspect the SingleCellExperiment object.

Preprocessing

Computation of QC metrics

Identify mitochondrial genes.

Calculate the Proportion of Mitochondrial Reads

Inspect the QC metrics

Filtering of low-quality cells

Visualise QC Metrics.

Filter Out Low Quality Cells

Visualise QC Metrics After Filtering

Further Preprocessing

Log-normalisation

Model gene variance

Dimensionality Reduction

Principal Component Analysis

Perform the PCA

Visualise the PCA Results

Visualise the PCA Results - Dimensional Loadings

Visualise the PCA Results - DimPlot

Visualise the PCA Results - Feature Plots

Visualise the PCA Results - Heatmaps

Determine the ‘dimensionality’ of the dataset

Make an Elbow Plot

Cluster the Cells

Run non-linear dimensional reduction (UMAP/tSNE)

Perform Dimensional Reduction with UMAP

Visualise UMAP

Visualise the UMAP

Other ways to visualise clusters

Visualise Expression across Clusters with a Violin Plot

Find Marker Genes

Find Positive Markers for Every Cluster Compared to the Rest

Use findMarkers to Compare All Clusters

Find Markers of Cluster 2

Use FindMarkers on a Single Cluster

Identify Cell Types

Make Violin Plots of B Cell Markers

Make Violin Plots of T Cell Markers

Colour the UMAP Plot by Canonical Marker Expression

Use Violin Plots to Compare Expression of Canonical Markers by Cluster

Rename Clusters with Cell Types

Revisualise the UMAP with Cell Type Annotations

Session information

Inspect the `SingleCellExperiment` object.

Use `findMarkers` to Compare All Clusters

Use `FindMarkers` on a Single Cluster