Introduction
Single-cell RNA sequencing (scRNA-seq) is a powerful and promising class of high-throughput assays that enable researchers to measure genome-wide transcription levels at the resolution of single cells. To properly account for features specific to scRNA-seq, such as zero inflation and high levels of technical noise, several novel statistical methods have been developed to tackle questions that include normalization, dimensionality reduction, clustering, the inference of cell lineages and pseudotimes, and the identification of differentially expressed (DE) genes. While each individual method is useful on its own for addressing a specific question, there is an increasing need for workflows that integrate these tools to yield a seamless scRNA-seq data analysis pipeline. This is all the more true, with novel sequencing technologies that allow an increasing number of cells to be sequenced in each run. For example, the Chromium Single Cell 3’ Solution was recently used to sequence and profile about 1.3 million cells from embryonic mouse brains.
scRNA-seq low-level analysis workflows have already been developed, with useful methods for quality control (QC), exploratory data analysis (EDA), pre-processing, normalization, and visualization. The workflow described in (Lun, McCarthy, and Marioni 2016) and the package scater (D. McCarthy et al. 2017) are such examples based on open-source R software packages from the Bioconductor Project (Huber et al. 2015). In these workflows, single-cell expression data are organized in objects of the SCESet class allowing integrated analysis. However, these workflows are mostly used to prepare the data for further downstream analysis and do not focus on steps such as cell clustering and lineage inference.
Here, we propose an integrated workflow for dowstream analysis, with the following four main steps: (1) dimensionality reduction accounting for zero inflation and over-dispersion and adjusting for gene and cell-level covariates, using the zinbwave Bioconductor package; (2) robust and stable cell clustering using resampling-based sequential ensemble clustering, as implemented in the clusterExperiment Bioconductor package; (3) inference of cell lineages and ordering of the cells by developmental progression along lineages, using the slingshot R package; and (4) DE analysis along lineages. Throughout the workflow, we use a single SummarizedExperiment object to store the scRNA-seq data along with any gene or cell-level metadata available from the experiment.
Analysis of olfactory stem cell differentiation using scRNA-seq data
Overview
This workflow is illustrated using data from a scRNA-seq study of stem cell differentiation in the mouse olfactory epithelium (OE) (Fletcher et al. 2017). The olfactory epithelium contains mature olfactory sensory neurons (mOSN) that are continuously renewed in the epithelium via neurogenesis through the differentiation of globose basal cells (GBC), which are the actively proliferating cells in the epithelium. When a severe injury to the entire tissue happens, the olfactory epithelium can regenerate from normally quiescent stem cells called horizontal basal cells (HBC), which become activated to differentiate and reconstitute all major cell types in the epithelium.
The scRNA-seq dataset we use as a case study was generated to study the differentitation of HBC stem cells into different cell types present in the olfactory epithelium. To map the developmental trajectories of the multiple cell lineages arising from HBCs, scRNA-seq was performed on FACS-purified cells using the Fluidigm C1 microfluidics cell capture platform followed by Illumina sequencing. The expression level of each gene in a given cell was quantified by counting the total number of reads mapping to it. Cells were then assigned to different lineages using a statistical analysis pipeline analogous to that in the present workflow. Finally, results were validated experimentally using in vivo lineage tracing. Details on data generation and statistical methods are available in (Fletcher et al. 2017; Risso et al. 2017; K. Street et al. 2017).
It was found that the first major bifurcation in the HBC lineage trajectory occurs prior to cell division, producing either mature sustentacular (mSUS) cells or GBCs. Then, the GBC lineage, in turn, branches off to give rise to mOSN, microvillous (MV) cells, and cells of the Bowman gland (Figure @ref(fig:stemcelldiff)). In this workflow, we describe a sequence of steps to recover the lineages found in the original study, starting from the genes x cells matrix of raw counts publicly-available at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE95601.
Package versions
The following packages are needed.
# Bioconductor
library(BiocParallel)
library(clusterExperiment)
library(scone)
library(zinbwave)
# GitHub
library(slingshot)
# CRAN
library(doParallel)
library(gam)
library(RColorBrewer)
set.seed(20)
Note that in order to successfully run the workflow, we need the devel versions of the Bioconductor packages scone (>=1.1.2), zinbwave (>=0.99.6), and clusterExperiment (>=1.3.2). We recommend running Bioconductor 3.6 (currently the devel version; see https://www.bioconductor.org/developers/how-to/useDevel/).
Parallel computing
For the workshop, we run the workflow in serial mode and do not run the time consuming functions zinbwave and RSEC. When running the workflow from scratch, we recommend running the workflow in parallel. See chunks below.
register(SerialParam())
NCORES <- 2
mysystem = Sys.info()[["sysname"]]
if (mysystem == "Darwin"){
registerDoParallel(NCORES)
register(DoparParam())
}else if (mysystem == "Linux"){
register(bpstart(MulticoreParam(workers=NCORES)))
}else{
print("Please change this to allow parallel computing on your computer.")
register(SerialParam())
}
Pre-processing
Counts for all genes in each cell were obtained from NCBI Gene Expression Omnibus (GEO), with accession number GSE95601. Before filtering, the dataset has 849 cells and 28,361 detected genes (i.e., genes with non-zero read counts).
data_dir <- "../data/"
if (!dir.exists(data_dir)) system(sprintf('mkdir %s', data_dir))
urls = c("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE95601&format=file&file=GSE95601%5FoeHBCdiff%5FCufflinks%5FeSet%2ERda%2Egz",
"https://raw.githubusercontent.com/rufletch/p63-HBC-diff/master/ref/oeHBCdiff_clusterLabels.txt")
if(!file.exists(paste0(data_dir, "GSE95601_oeHBCdiff_Cufflinks_eSet.Rda"))) {
if (!dir.exists(data_dir)) system(sprintf('mkdir %s', data_dir))
download.file(urls[1], paste0(data_dir, "GSE95601_oeHBCdiff_Cufflinks_eSet.Rda.gz"))
R.utils::gunzip(paste0(data_dir, "GSE95601_oeHBCdiff_Cufflinks_eSet.Rda.gz"))
assayData(Cufflinks_eSet)$exprs = NULL
assayData(Cufflinks_eSet)$fpkm_table = NULL
assayData(Cufflinks_eSet)$tpm_table = NULL
save(Cufflinks_eSet, file='data/GSE95601_oeHBCdiff_Cufflinks_eSet_reduced.Rda')
}
if(!file.exists(paste0(data_dir, "oeHBCdiff_clusterLabels.txt"))) {
download.file(urls[2], paste0(data_dir, "oeHBCdiff_clusterLabels.txt"))
}
load(paste0(data_dir, "GSE95601_oeHBCdiff_Cufflinks_eSet_reduced.Rda"))
# Count matrix
E <- assayData(Cufflinks_eSet)$counts_table
# Remove undetected genes
E <- na.omit(E)
E <- E[rowSums(E)>0,]
dim(E)
## [1] 28361 849
We remove the ERCC spike-in sequences and the CreER gene, as the latter corresponds to the estrogen receptor fused to Cre recombinase (Cre-ER), which is used to activate HBCs into differentiation following injection of tamoxifen (see (Fletcher et al. 2017) for details).
# Remove ERCC and CreER genes
cre <- E["CreER",]
ercc <- E[grep("^ERCC-", rownames(E)),]
E <- E[grep("^ERCC-", rownames(E), invert = TRUE), ]
E <- E[-which(rownames(E)=="CreER"), ]
dim(E)
## [1] 28284 849
Throughout the workflow, we use the class SummarizedExperiment to keep track of the counts and their associated metadata within a single object. The cell-level metadata contain quality control measures, sequencing batch ID, and cluster and lineage labels from the original publication (Fletcher et al. 2017). Cells with a cluster label of -2 were not assigned to any cluster in the original publication.
# Extract QC metrics
qc <- as.matrix(protocolData(Cufflinks_eSet)@data)[,c(1:5, 10:18)]
qc <- cbind(qc, CreER = cre, ERCC_reads = colSums(ercc))
# Extract metadata
batch <- droplevels(pData(Cufflinks_eSet)$MD_c1_run_id)
bio <- droplevels(pData(Cufflinks_eSet)$MD_expt_condition)
clusterLabels <- read.table(paste0(data_dir, "oeHBCdiff_clusterLabels.txt"),
sep = "\t", stringsAsFactors = FALSE)
m <- match(colnames(E), clusterLabels[, 1])
# Create metadata data.frame
metadata <- data.frame("Experiment" = bio,
"Batch" = batch,
"publishedClusters" = clusterLabels[m,2],
qc)
# Symbol for cells not assigned to a lineage in original data
metadata$publishedClusters[is.na(metadata$publishedClusters)] <- -2
se <- SummarizedExperiment(assays = list(counts = E),
colData = metadata)
se
## class: SummarizedExperiment
## dim: 28284 849
## metadata(0):
## assays(1): counts
## rownames(28284): Xkr4 LOC102640625 ... Ggcx.1 eGFP
## rowData names(0):
## colnames(849): OEP01_N706_S501 OEP01_N701_S501 ... OEL23_N704_S503
## OEL23_N703_S502
## colData names(19): Experiment Batch ... CreER ERCC_reads
Using the Bioconductor R package scone, we remove low-quality cells according to the quality control filter implemented in the function metric_sample_filter and based on the following criteria (Figure @ref(fig:scone)): (1) Filter out samples with low total number of reads or low alignment percentage and (2) filter out samples with a low detection rate for housekeeping genes. See the scone vignette for details on the filtering procedure.
# QC-metric-based sample-filtering
data("housekeeping")
hk = rownames(se)[toupper(rownames(se)) %in% housekeeping$V1]
mfilt <- metric_sample_filter(assay(se),
nreads = colData(se)$NREADS,
ralign = colData(se)$RALIGN,
pos_controls = rownames(se) %in% hk,
zcut = 3, mixture = FALSE,
plot = TRUE)
# Simplify to a single logical
mfilt <- !apply(simplify2array(mfilt[!is.na(mfilt)]), 1, any)
se <- se[, mfilt]
dim(se)
## [1] 28284 747
After sample filtering, we are left with 747 good quality cells.
Finally, for computational efficiency, we retain only the 1,000 most variable genes. This seems to be a reasonnable choice for the illustrative purpose of this workflow, as we are able to recover the biological signal found in the published analysis ((Fletcher et al. 2017)). In general, however, we recommend care in selecting a gene filtering scheme, as an appropriate choice is dataset-dependent.
# Filtering to top 1,000 most variable genes
vars <- rowVars(log1p(assay(se)))
names(vars) <- rownames(se)
vars <- sort(vars, decreasing = TRUE)
core <- se[names(vars)[1:1000],]
Dataset structure
Overall, after the above pre-processing steps, our dataset has 1,000 genes and 747 cells.
core
## class: SummarizedExperiment
## dim: 1000 747
## metadata(0):
## assays(1): counts
## rownames(1000): Cbr2 Cyp2f2 ... Rnf13 Atp7b
## rowData names(0):
## colnames(747): OEP01_N706_S501 OEP01_N701_S501 ... OEL23_N704_S503
## OEL23_N703_S502
## colData names(19): Experiment Batch ... CreER ERCC_reads
Metadata for the cells are stored in the slot colData from the SummarizedExperiment object. Cells were processed in 18 different batches.
batch <- colData(core)$Batch
col_batch = c(brewer.pal(9, "Set1"), brewer.pal(8, "Dark2"),
brewer.pal(8, "Accent")[1])
names(col_batch) = unique(batch)
table(batch)
## batch
## GBC08A GBC08B GBC09A GBC09B P01 P02 P03A P03B P04 P05
## 39 40 35 22 31 48 51 40 20 23
## P06 P10 P11 P12 P13 P14 Y01 Y04
## 51 40 50 50 60 47 58 42
In the original work (Fletcher et al. 2017), cells were clustered into 14 different clusters, with 151 cells not assigned to any cluster (i.e., cluster label of -2).
publishedClusters <- colData(core)[, "publishedClusters"]
col_clus <- c("transparent", "#1B9E77", "antiquewhite2", "cyan", "#E7298A",
"#A6CEE3", "#666666", "#E6AB02", "#FFED6F", "darkorchid2",
"#B3DE69", "#FF7F00", "#A6761D", "#1F78B4")
names(col_clus) <- sort(unique(publishedClusters))
table(publishedClusters)
## publishedClusters
## -2 1 2 3 4 5 7 8 9 10 11 12 14 15
## 151 90 25 54 35 93 58 27 74 26 21 35 26 32
Note that there is partial nesting of batches within clusters (i.e., cell type), which could be problematic when correcting for batch effects in the dimensionality reduction step below.
table(data.frame(batch = as.vector(batch),
cluster = publishedClusters))
## cluster
## batch -2 1 2 3 4 5 7 8 9 10 11 12 14 15
## GBC08A 3 0 2 12 9 0 0 0 0 0 2 0 2 9
## GBC08B 8 0 7 5 3 0 0 0 1 2 3 0 5 6
## GBC09A 6 0 1 5 8 0 0 0 1 1 0 0 6 7
## GBC09B 12 0 2 1 3 0 0 0 1 0 0 0 3 0
## P01 7 0 2 4 3 15 0 0 0 0 0 0 0 0
## P02 5 2 0 9 3 15 3 3 2 3 0 2 1 0
## P03A 15 3 0 2 0 12 2 9 4 2 0 2 0 0
## P03B 9 1 2 1 1 11 1 2 8 1 1 2 0 0
## P04 8 0 0 0 0 9 1 0 1 1 0 0 0 0
## P05 3 0 0 0 1 11 3 0 1 0 2 2 0 0
## P06 12 1 2 3 0 8 2 4 8 4 1 2 2 2
## P10 7 3 1 4 0 3 5 8 1 0 2 5 0 1
## P11 6 2 1 1 0 1 5 1 22 3 1 6 0 1
## P12 10 0 2 0 0 4 10 0 8 2 3 6 4 1
## P13 13 1 2 4 0 4 15 0 4 5 6 1 3 2
## P14 9 0 0 1 2 0 11 0 12 2 0 7 0 3
## Y01 8 46 1 1 2 0 0 0 0 0 0 0 0 0
## Y04 10 31 0 1 0 0 0 0 0 0 0 0 0 0
Normalization and dimensionality reduction: ZINB-WaVE
In scRNA-seq analysis, dimensionality reduction is often used as a preliminary step prior to downstream analyses, such as clustering, cell lineage and pseudotime ordering, and the identification of DE genes. This allows the data to become more tractable, both from a statistical (cf. curse of dimensionality) and computational point of view. Additionally, technical noise can be reduced while preserving the often intrinsically low-dimensional signal of interest (Dijk et al. 2017; Pierson and Yau 2015; Risso et al. 2017).
Here, we perform dimensionality reduction using the zero-inflated negative binomial-based wanted variation extraction (ZINB-WaVE) method implemented in the Bioconductor R package zinbwave. The method fits a ZINB model that accounts for zero inflation (dropouts), over-dispersion, and the count nature of the data. The model can include a cell-level intercept, which serves as a global-scaling normalization factor. The user can also specify both gene-level and cell-level covariates. The inclusion of observed and unobserved cell-level covariates enables normalization for complex, non-linear effects (often referred to as batch effects), while gene-level covariates may be used to adjust for sequence composition effects (e.g., gene length and GC-content effects). A schematic view of the ZINB-WaVE model is provided in Figure @ref(fig:zinbschema). For greater detail about the ZINB-WaVE model and estimation procedure, please refer to the original manuscript (Risso et al. 2017).
As with most dimensionality reduction methods, the user needs to specify the number of dimensions for the new low-dimensional space. Here, we use K = 50 dimensions and adjust for batch effects via the matrix X. Note that if the users include more genes in the analysis, it may be preferable to reduce K to achieve a similar computational time.
print(system.time(se <- zinbwave(core, K = 50, X = "~ Batch",
residuals = TRUE,
normalizedValues = TRUE)))
save(se, file = 'se_after_zinbwave.rda')
load(sprintf('%sse_after_zinbwave.rda', data_dir))
Normalization
The function zinbwave returns a SummarizedExperiment object that includes normalized expression measures, defined as deviance residuals from the fit of the ZINB-WaVE model with user-specified gene- and cell-level covariates. Such residuals can be used for visualization purposes (e.g., in heatmaps, boxplots). Note that, in this case, the low-dimensional matrix W is not included in the computation of residuals to avoid the removal of the biological signal of interest.
norm <- assays(se)$normalizedValues
norm[1:3,1:3]
## OEP01_N706_S501 OEP01_N701_S501 OEP01_N707_S507
## Cbr2 4.557371 4.375069 -4.142697
## Cyp2f2 4.321644 4.283266 4.090283
## Gstm1 4.796498 4.663366 4.416324
As expected, the normalized values no longer exhibit batch effects (Figure @ref(fig:boxplotNorm)).
norm_order <- norm[, order(as.numeric(batch))]
col_order <- col_batch[batch[order(as.numeric(batch))]]
boxplot(norm_order, col = col_order, staplewex = 0, outline = 0,
border = col_order, xaxt = "n", ylab="Expression measure")
abline(h=0)
The principal component analysis (PCA) of the normalized values shows that, as expected, cells do not cluster by batch but by the original clusters (Figure @ref(fig:pcanorm)). Overall, it seems that normalization was effective at removing batch effects without removing biological signal, in spite of the partial nesting of batches within clusters.
pca <- prcomp(t(norm))
par(mfrow = c(1,2))
plot(pca$x, col = col_batch[batch], pch = 20, main = "")
plot(pca$x, col = col_clus[as.character(publishedClusters)], pch = 20, main = "")
