2.1 Introducing the calculateQCMetrics function

The calculateQCMetrics function computes a number of quality control metrics for each cell and feature, stored in the colData and rowData respectively. By default, the QC metrics are computed from the count data, but this can be changed through the exprs_values argument.

example_sce <- calculateQCMetrics(example_sce)
colnames(colData(example_sce))

##  [1] "Cell"                           "Mutation_Status"               
##  [3] "Cell_Cycle"                     "Treatment"                     
##  [5] "is_cell_control"                "total_features_by_counts"      
##  [7] "log10_total_features_by_counts" "total_counts"                  
##  [9] "log10_total_counts"             "pct_counts_in_top_50_features" 
## [11] "pct_counts_in_top_100_features" "pct_counts_in_top_200_features"
## [13] "pct_counts_in_top_500_features" "total_features"                
## [15] "log10_total_features"           "pct_counts_top_50_features"    
## [17] "pct_counts_top_100_features"    "pct_counts_top_200_features"   
## [19] "pct_counts_top_500_features"

colnames(rowData(example_sce))

## [1] "is_feature_control"    "mean_counts"           "log10_mean_counts"    
## [4] "n_cells_by_counts"     "pct_dropout_by_counts" "total_counts"         
## [7] "log10_total_counts"    "n_cells_counts"        "pct_dropout_counts"

Control sets can be defined for features (e.g., spike-in transcripts, mitochondrial genes) or cells (e.g., empty wells, visually damaged cells). The function will subsequently compute metrics for each control set, e.g., the proportion of counts assigned to spike-in transcripts.

example_sce <- calculateQCMetrics(example_sce, 
    feature_controls = list(ERCC = 1:20, mito = 500:1000),
    cell_controls = list(empty = 1:5, damaged = 31:40))

all_col_qc <- colnames(colData(example_sce))
all_col_qc <- all_col_qc[grep("ERCC", all_col_qc)]

Metrics are also computed for the union of all control sets, and for the set of all features/cells that are not marked as controls. These sets are labelled as "feature_control" and "endogenous" for features, and "cell_control" and "non_control" for cells. Users should avoid using these names for their own sets.

2.2 Cell-level QC metrics

We refer users to the ?calculateQCMetrics for a full list of the computed cell-level metrics. However, we will describe some of the more popular ones here:

• total_counts: total number of counts for the cell (i.e., the library size).
• total_features_by_counts: the number of features for the cell that have counts above the detection limit (default of zero).
• pct_counts_X: percentage of all counts that come from the feature control set named X.

If exprs_values is set to something other than "counts", the names of the metrics will be changed by swapping "counts" for whatever named assay was used.

2.3 Feature-level QC metrics

Feature-level metrics include:

• mean_counts: the mean count of the gene/feature.
• pct_dropout_by_counts: the percentage of cells with counts of zero for each gene.
• pct_counts_Y: percentage of all counts that come from the cell control set named Y.

Again, if a different exprs_values was used, the names of the metrics will change correspondingly.

3.1 Examining the most expressed features

We look at a plot that shows the top 50 (by default) most-expressed features. Each row in the plot below corresponds to a gene, and each bar corresponds to the expression of a gene in a single cell. The circle indicates the median expression of each gene, with which genes are sorted. By default, “expression” is defined using the feature counts (if available), but other expression values can be used instead by changing exprs_values.

plotQC(example_sce, type = "highest-expression", exprs_values = "counts")

We expect to see the “usual suspects”, i.e., mitochondrial genes, actin, ribosomal protein, MALAT1. A few spike-in transcripts may also be present here, though if all of the spike-ins are in the top 50, it suggests that too much spike-in RNA was added. A large number of pseudo-genes or predicted genes may indicate problems with alignment.

3.2 Frequency of expression as a function of the mean

Another useful plot is that of the frequency of expression (i.e., number of cells with non-zero expression) against the mean. These two metrics should be positively correlated with each other for most genes.

plotQC(example_sce, type = "exprs-freq-vs-mean")

Outliers from the trend may warrant further investigation. For example, alignment errors for pseudo-genes of highly-expressed genes will result in features with low means that are expressed in all cells. Conversely, PCR amplification biases (or the presence of rare populations) may result in genes with very high means that are expressed in very few cells.

3.3 Percentage of counts assigned to feature controls

A particularly useful plot for cell-level QC involves percentage of expression in feature controls against the total number of expressed features. These two metadata variables can be plotted against each other as shown below. We take advantage of the ggplot2 semantics to fine-tune the plot aesthetics and to add a smoothing curve:

plotColData(example_sce, x = "total_features_by_counts",
    y = "pct_counts_feature_control", colour = "Mutation_Status") +
    theme(legend.position = "top") +
    stat_smooth(method = "lm", se = FALSE, size = 2, fullrange = TRUE)

Well-behaved cells should have a large number of expressed features and and a low percentage of expression from feature controls. High percentage expression from feature controls and few expressed features are indicative of blank and failed cells.

3.4 Cumulative expression plot

The plotScater method plots the cumulative proportion of each cell’s library assigned to the top highest-expressed features (default 500). This type of plot visualizes differences in expression distributions for different cells, in the same manner as per-sample boxplots for microarray or bulk RNA-seq data. Cumulative expression plots are more effective for single-cell data where it is not easy to examine hundreds or thousands of boxplots at once.

With this function, we can split up the cells based on colData variables to examine differences in the expression distributions between cells. By default, the plot method will try to use count values for the plot, but other data available in the assays slot of the object can be used by specifying exprs_values.

plotScater(example_sce, block1 = "Mutation_Status", block2 = "Treatment",
     colour_by = "Cell_Cycle", nfeatures = 300, exprs_values = "counts")

This approach can allow users to identify large differences in expression distributions across different experimental blocks (e.g., processing batches).

3.5 Plate position plot

For plate-based experiments, it is useful to see how expression or factors vary with the position of cell on the plate. This can be visualized using the plotPlatePosition function:

example_sce2 <- example_sce
example_sce2$plate_position <- paste0(
     rep(LETTERS[1:5], each = 8), 
     rep(formatC(1:8, width = 2, flag = "0"), 5)
)
plotPlatePosition(example_sce2, colour_by = "Gene_0001",
    by_exprs_values = "counts") 

Systematic trends in expression with the plate position may indicate that there were issues with processing. The same approach can be used with experimental factors to determine whether cells are appropriately randomized across the plate.

3.6 Other quality control plots

Two feature metadata variables can be easily plotted against each other using the plotFeatureData function:

plotRowData(example_sce, x = "n_cells_by_counts", y = "mean_counts")

The multiplot function also allows multiple plots to be generated on the same page, as demonstrated below.

p1 <- plotColData(example_sce, x = "total_counts", 
    y = "total_features_by_counts")
p2 <- plotColData(example_sce, x = "pct_counts_feature_control",
    y = "total_features_by_counts")
p3 <- plotColData(example_sce, x = "pct_counts_feature_control",
    y = "pct_counts_top_50_features")
multiplot(p1, p2, p3, cols = 3)

This is especially useful for side-by-side comparisons between control sets, as demonstrated below for the plot of highest-expressing features. A plot for non-control cells is shown on the left while the plot for the controls is shown on the right.

p1 <- plotQC(example_sce[, !example_sce$is_cell_control],
    type = "highest-expression")
p2 <- plotQC(example_sce[, example_sce$is_cell_control],
    type = "highest-expression")
multiplot(p1, p2, cols = 2)

4.1 By cells

example_sce <- example_sce[,1:40]

filter(example_sce, Treatment == "treat1")

## class: SingleCellExperiment 
## dim: 2000 27 
## metadata(0):
## assays(1): counts
## rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
## rowData names(47): is_feature_control is_feature_control_ERCC ...
##   n_cells_counts_damaged pct_dropout_counts_damaged
## colnames(27): Cell_001 Cell_002 ... Cell_037 Cell_039
## colData names(73): Cell Mutation_Status ...
##   pct_counts_top_200_features_mito pct_counts_top_500_features_mito
## reducedDimNames(0):
## spikeNames(0):

keep.total <- example_sce$total_counts > 1e5
keep.n <- example_sce$total_features_by_counts > 500
filtered <- example_sce[,keep.total & keep.n]
dim(filtered)

## [1] 2000   37

keep.total <- isOutlier(example_sce$total_counts, nmads=3, 
    type="lower", log=TRUE)
filtered <- example_sce[,keep.total]

example_sce <- runPCA(example_sce, use_coldata = TRUE,
    detect_outliers = TRUE)
plotReducedDim(example_sce, use_dimred="PCA_coldata")

summary(example_sce$outlier)

##    Mode   FALSE 
## logical      40

4.2 By features

It is common to filter out low-abundance features prior to further analyses. This is easily achieved by row subsetting of the SingleCellExperiment object. In the example below, genes are only retained if they are expressed in four or more cells:

keep_feature <- nexprs(example_sce, byrow=TRUE) >= 4
example_sce <- example_sce[keep_feature,]
dim(example_sce)

## [1] 1753   40

Other filtering can be done using existing annotation. For example, ribosomal protein genes and predicted genes can be identified (and removed) using regular expressions or biotype information. Such genes are often (but not always) uninteresting for characterizing population heterogeneity.

6.1 Scaling normalization

Scaling normalization accounts for cell-specific biases that scale expression up or down for all genes in a particular cell, e.g., coverage or capture efficiency. The simplest approach to scaling normalization defines the size factors from the scaled library sizes of all cells. is done so that the mean size factor is equal to unity, ensuring that the normalized values are on the same scale as the original counts.

sizeFactors(example_sce) <- librarySizeFactors(example_sce)
summary(sizeFactors(example_sce))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1463  0.6609  0.8112  1.0000  1.2533  2.7356

Log-transformed normalized expression values can then be computed with normalize, which stores the output in the "logcounts" slot.

example_sce <- normalize(example_sce)

While simple, library size normalization does not account for composition biases that are often present in high-throughput sequencing data. It also fails to account for differences in the biases affecting spike-in transcripts. We strongly suggest using the computeSumFactors and computeSpikeFactors functions from the scran package.

6.2 Batch correction

Batch correction accounts for systematic differences in expression between cells in different batches. Unlike scaling biases, these are usually constant across all cells in a given batch but different for each gene.

Batch effects can be regressed out by using the removeBatchEffect function from the limma package. This applies a linear model, usually on the log-expression values to avoid issues with the mean-variance relationship. To illustrate:

library(limma)
batch <- rep(1:2, each=20)
corrected <- removeBatchEffect(logcounts(example_sce), block=batch)
assay(example_sce, "corrected_logcounts") <- corrected

Factors of interest can be included in design to avoid regressing them out. This is necessary when they are not orthogonal to the block. However, this assumes that your model is fully specified, which may not be possible when the factors of interest are unknown. In such cases, an alternative method is to use the mnnCorrect approach from scran.