Single Cells with edgeR

This vignette demonstrates usage of Glimma on a single cell dataset. The data here comes from brain cells from the Zeisel A et al. (2015) study on mouse brain single cells. We will use the MDS plot to perform unsupervised clustering of the cells. A pseudo-bulk single cell aggregation approach with edgeR will be used to test for differential expression, and two styles of MA plots will be used to investigate the results.

This is a simplified workflow not intended to represent best practice, but to produce reasonable looking plots in the minimal amount of code. Please refer to a resource such Orchestrating Single-Cell Analysis with Bioconductor (OSCA) for appropriate workflows for analysis.

We start by loading in the data using the scRNAseq package.

library(scRNAseq)
library(scater)
library(scran)
library(Glimma)
library(edgeR)

sce <- ZeiselBrainData(ensembl=TRUE)
#> snapshotDate(): 2022-04-19
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> snapshotDate(): 2022-04-19
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> snapshotDate(): 2022-04-21
#> loading from cache
#> require("ensembldb")
#> Warning: Unable to map 1565 of 20006 requested IDs.

Once the data is loaded in we follow the OSCA procedure for identifying highly variable genes for creating a multi-dimensional scaling (MDS) plot. We use the functions provided by scran to identify the most highly variable genes rather than the algorithm within glimmaMDS, as scran is tailored towards single cells.

sce <- logNormCounts(sce)

var_mod <- modelGeneVar(sce)
hvg_genes <- getTopHVGs(var_mod, n=500)
hvg_sce <- sce[hvg_genes, ]

hvg_sce <- logNormCounts(hvg_sce)

Choosing to colour the MDS plot using level1class reveals separation between cell types.

glimmaMDS(
    exprs(hvg_sce),
    groups = colData(hvg_sce)
)