scRNA-seq workflow: Macosko et al.

1 Import data

The scRNAseq package provides convenient access to several datasets. See the package Bioconductor page for more information.

# install BiocManager package if not installed yet.
# BiocManager is the package installer for Bioconductor software.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# install packages if not yet installed.
pkgs <- c("SingleCellExperiment", "DropletUtils", "scRNAseq", "scater", "scuttle", "scran", "BiocSingular", "scDblFinder", "glmpca", "uwot")
notInstalled <- pkgs[!pkgs %in% installed.packages()[,1]]
if(length(notInstalled) > 0){
  BiocManager::install(notInstalled)
}

# Code below might ask you to create an ExperimentHub directory. 
# Type 'yes' and hit Enter, to allow this.
suppressPackageStartupMessages(library(scRNAseq))
sce <- MacoskoRetinaData()

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

## see ?scRNAseq and browseVignettes('scRNAseq') for documentation

## loading from cache

2 A SingleCellExperiment object

sce

## class: SingleCellExperiment 
## dim: 24658 49300 
## metadata(0):
## assays(1): counts
## rownames(24658): KITL TMTC3 ... 1110059M19RIK GM20861
## rowData names(0):
## colnames(49300): r1_GGCCGCAGTCCG r1_CTTGTGCGGGAA ... p1_TAACGCGCTCCT
##   p1_ATTCTTGTTCTT
## colData names(2): cell.id cluster
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

2.1 Accessing data from a SingleCellExperiment object

Please see Figure 4.1 in OSCA for an overview of a SingleCellExperiment object.

# Data: assays
assays(sce)

## List of length 1
## names(1): counts

assays(sce)$counts[1:5, 1:5]

## 5 x 5 sparse Matrix of class "dgCMatrix"
##               r1_GGCCGCAGTCCG r1_CTTGTGCGGGAA r1_GCGCAACTGCTC r1_GATTGGGAGGCA
## KITL                        .               .               1               .
## TMTC3                       3               .               .               .
## CEP290                      1               3               .               2
## 4930430F08RIK               2               1               2               .
## 1700017N19RIK               .               .               .               .
##               r1_CCTCCTAGTTGG
## KITL                        .
## TMTC3                       2
## CEP290                      1
## 4930430F08RIK               1
## 1700017N19RIK               .

# Feature metadata: rowData
rowData(sce) # empty for now

## DataFrame with 24658 rows and 0 columns

# Cell metadata: colData
colData(sce)

## DataFrame with 49300 rows and 2 columns
##                         cell.id   cluster
##                     <character> <integer>
## r1_GGCCGCAGTCCG r1_GGCCGCAGTCCG         2
## r1_CTTGTGCGGGAA r1_CTTGTGCGGGAA         2
## r1_GCGCAACTGCTC r1_GCGCAACTGCTC         2
## r1_GATTGGGAGGCA r1_GATTGGGAGGCA         2
## r1_CCTCCTAGTTGG r1_CCTCCTAGTTGG        NA
## ...                         ...       ...
## p1_TCAAAAGCCGGG p1_TCAAAAGCCGGG        24
## p1_ATTAAGTTCCAA p1_ATTAAGTTCCAA        34
## p1_CTGTCTGAGACC p1_CTGTCTGAGACC         2
## p1_TAACGCGCTCCT p1_TAACGCGCTCCT        24
## p1_ATTCTTGTTCTT p1_ATTCTTGTTCTT        24

# Reduced dimensions: reducedDims
reducedDims(sce) # empty for now

## List of length 0
## names(0):

2.2 Creating a new SingleCellExperiment object

sceNew <- SingleCellExperiment(assays = list(counts = assays(sce)$counts))
sceNew

## class: SingleCellExperiment 
## dim: 24658 49300 
## metadata(0):
## assays(1): counts
## rownames(24658): KITL TMTC3 ... 1110059M19RIK GM20861
## rowData names(0):
## colnames(49300): r1_GGCCGCAGTCCG r1_CTTGTGCGGGAA ... p1_TAACGCGCTCCT
##   p1_ATTCTTGTTCTT
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

rm(sceNew)

2.3 Storing (meta)data in a SingleCellExperiment object

fakeGeneNames <- paste0("gene", 1:nrow(sce))
rowData(sce)$fakeName <- fakeGeneNames
head(rowData(sce))

## DataFrame with 6 rows and 1 column
##                  fakeName
##               <character>
## KITL                gene1
## TMTC3               gene2
## CEP290              gene3
## 4930430F08RIK       gene4
## 1700017N19RIK       gene5
## MGAT4C              gene6

# Remove again by setting to NULL
rowData(sce)$fakeName <- NULL

assays(sce)$logCounts <- log1p(assays(sce)$counts)
assays(sce)

## List of length 2
## names(2): counts logCounts

assays(sce)$logCounts[1:5, 1:5]

## 5 x 5 sparse Matrix of class "dgCMatrix"
##               r1_GGCCGCAGTCCG r1_CTTGTGCGGGAA r1_GCGCAACTGCTC r1_GATTGGGAGGCA
## KITL                .               .               0.6931472        .       
## TMTC3               1.3862944       .               .                .       
## CEP290              0.6931472       1.3862944       .                1.098612
## 4930430F08RIK       1.0986123       0.6931472       1.0986123        .       
## 1700017N19RIK       .               .               .                .       
##               r1_CCTCCTAGTTGG
## KITL                .        
## TMTC3               1.0986123
## CEP290              0.6931472
## 4930430F08RIK       0.6931472
## 1700017N19RIK       .

assays(sce)$logCounts <- NULL

3 Filtering non-informative genes

4 Quality control

4.1 Calculate QC variables

Note that the character string “^MT-” in the rownames of the sce object indicates mitochondrial genes.

4.2 EDA

High-quality cells should have many features expressed, and a low contribution of mitochondrial genes. Here, we see that several cells have a very low number of expressed genes, and where most of the molecules are derived from mitochondrial genes. This indicates likely damaged cells, presumably because of loss of cytoplasmic RNA from perforated cells, so we’d want to remove these for the downstream analysis.

4.3 QC using adaptive thresholds

Below, we remove cells that are outlying with respect to

1. A low sequencing depth (number of UMIs);
2. A low number of genes detected;
3. A high percentage of reads from mitochondrial genes.

Also make plots to assess the impact of the removal of .

4.4 Identifying and removing empty droplets

Note that the removal of cells with low sequencing depth using the adaptive threshold procedure above is a way of removing empty droplets. Other approaches are possible, e.g., removing cells by statistical testing using emtpyDrops. This does require us to specify a lower bound on the total number of UMIs, below which all cells are considered to correspond to empty droplets. This lower bound may not be trivial to derive, but the barcodeRanks function can be useful to identify an elbow/knee point.

4.5 Identifying and removing doublets

We will use scDblFinder to detect doublet cells.

5 Normalization

For normalization, the size factors \(s_i\) computed here are simply scaled library sizes: \[ N_i = \sum_g Y_{gi} \] \[ s_i = N_i / \bar{N}_i \]

6 Feature selection

7 Dimension reduction

Note that, below, we color the cells using the known, true cell type label as defined in the metadata, to empirically evaluate the dimension reduction. In reality, we don’t know this yet at this stage.

7.1 Linear dimension reduction: PCA

We are able to recover quite some structure. However, many cell populations remain obscure, and the plot is overcrowded.

7.1.1 PCA without feature selection

7.2 A generalization of PCA for exponential family distributions.

7.3 Non-linear dimension reduction: UMAP

8 Clustering

# Build a shared nearest-neighbor graph from PCA space
#g <- buildSNNGraph(sce, use.dimred = 'PCA')

# Louvain clustering on the SNN graph, and add to sce
#colData(sce)$label <- factor(igraph::cluster_louvain(g)$membership)

# Visualization.
#plotUMAP(sce, colour_by="label")