3.2 SCESet: a data structure for single-cell expression data

The scater package is built around the SCESet class (Supplementary Fig. S1) which provides a sophisticated container for scRNA-seq data. This class inherits from the ExpressionSet class in Bioconductor’s Biobase package (Huber et al., 2015), which allows assay data (and multiple transformations thereof), gene or transcript metadata and sample metadata to be combined in a single object to empower robust analyses. While the ExpressionSet class is the basis of many microarray and bulk RNA-seq analysis methods in Bioconductor, extensions to the class design are necessary to support scRNA-seq data analyses. Specifically, the SCESet class adds slots to store a reduced-dimension representation of the expression profiles, to easily visualize the relationships between cells; cell–cell and gene–gene pairwise distance matrices, for clustering or regulatory network reconstruction; bootstrapped expression results (such as from kallisto), to gauge the accuracy of expression quantification; consensus clustering results, where cluster assignments for each cell are combined from different methods to improve reliability; information about feature controls (such as ERCC spike-ins), which is required in downstream steps such as normalization, QC and detection of highly variable genes; and several more (Supplementary Fig. S1). With these extra slots, SCESet objects can support analyses of scRNA-seq data that ExpressionSet cannot. In addition, extra data types such as FACS marker expression or epigenetic information can be easily stored in each SCESet object for integration with the single-cell expression profiles.

An SCESet data object can be easily subsetted by row or column to remove unwanted genes or cells, respectively, from all data and metadata fields stored in the object. Furthermore, data and metadata in multiple SCESet objects can be easily combined e.g. to incorporate cells from different experimental batches. SCESet objects can also be converted to other R data structures, or saved to disk in structured, shareable formats. Further details on the class, including its motivation and execution, are available in the Supplementary Case Study and the package documentation. All methods available in scater are applicable to instances of the SCESet class and exploit the availability and richness of (meta)data stored in each SCESet object.

3.3 Data pre-processing

An important initial step in scRNA-seq data processing is to quantify the expression level of genomic features such as transcripts or genes from the raw sequencing data. Approaches to expression quantification from raw reads are, in principle, the same for scRNA-seq as they are for bulk RNA-seq (Kanitz et al., 2015; Teng et al., 2016). Read counts obtained from conventional quantification methods such as HTSeq (Anders et al., 2015) and featureCounts (Liao et al., 2014) can be readily stored in an SCESet object and used in a scater workflow (Fig. 1). Another option is to use computationally efficient pseudoalignment methods such as kallisto and Salmon. This is especially appealing for large scRNA-seq datasets containing hundreds to tens of thousands of cells. To this end, scater also provides wrapper functions for kallisto and Salmon so that fast quantification of transcript-level expression can be managed completely within an R programming environment. A common subsequent step for these methods is to collapse transcript-level expression to gene-level expression. Exploiting the biomaRt R/Bioconductor package, scater provides a convenient function for using Ensembl annotations (Yates et al., 2016) to obtain gene-level expression values and gene or transcript annotations.

3.4 Data quality control

The scater package provides methods to compute relevant QC metrics for an SCESet object. Given a set of control genes and/or cells, a variety of QC metrics will be computed and returned to the object in a single call to the calculateQCMetrics function (see package documentation). Cell-specific QC metrics include the total count across all genes, the total number of expressed genes, and the percentage of counts allocated to control genes like spike-in transcripts or mitochondrial genes. These metrics are useful for identifying low-quality cells—for example, a high percentage of counts mapping to spike-ins typically indicates that a small amount of RNA was captured for the cell, suggesting protocol failure or death of the cell in processing that renders it unsuitable for downstream analyses. For each gene, QC metrics such as the average expression level and the proportion of cells in which the gene is expressed are computed. This can be used to identify low-abundance genes or genes with high dropout rates that should be filtered out prior to downstream analyses. All of these metrics are used by scater to construct QC plots to diagnose potential issues with data quality. This facilitates quality control which—despite attempts at automation (Ilicic et al., 2016)—still requires manual intervention to account for aspects of the data specific to each study. The package documentation provides full details of the QC metrics produced.

In scater, the default plot method for an SCESet object produces a cumulative expression plot (Fig. 2a). This plot describes how reads are distributed across genes, distinguishing between low-complexity libraries (where very few genes contain most of the counts) and their high-complexity counterparts (where counts are distributed more evenly across genes). For example, there is substantial variability in library complexity among cells in the case study dataset in Figure 2a. Some cells have profiles similar to the blank wells, suggesting that library preparation or sequencing failed for these cells and that the corresponding libraries should be removed prior to further analysis. Cell phenotype variables can be incorporated into these plots to highlight differences in expression distributions for different types of cells. For example, the curve for each cell is coloured by the type of well that produced the library (Fig. 2a), while cells can also be split into separate facets by library type to show more metadata variables simultaneously (see Supplementary Case Study). Cumulative expression plots should be favoured over boxplots as the default method for visualizing expression distributions across cells in a dataset, as the latter performs poorly at handling the long tail of low- and zero-expression observations in scRNA-seq data.

Open in a separate window

Fig. 2.

Different types of QC plots that can be generated with scater. (a) Cumulative expression plot showing the proportion of the library accounted for by the top 1–500 most highly expressed features. (b) PCA plot produced using a subset of the QC metrics computed with scater’s calculateQCMetrics function. (c) Plot of frequency of expression (percentage of cells in which the feature is deemed expressed) against mean expression level across cells. The vertical dotted line shows the median of the gene mean expression levels, and the horizontal dotted line indicates 50% frequency of expression. (d) Plot of the 20 most highly expressed features (computed according to the highest total read counts) across all cells in the dataset. For each feature, the circle represents the percentage of counts across all cells that correspond to that feature. The features are ordered by this value. The bars for each feature show the percentage of counts corresponding to the feature in each individual cell, providing a visualization of the distribution across cells. (e) Density plot showing the percentage of variance explained by a set of explanatory variables across all genes. Each individual plot is produced by a single call with either the function plot (a), plotPCA (b) or plotQC (c–e)

The plotPCA function implements an approach to automatic outlier detection using multivariate normal methods applied to the cell-level QC metrics (Ilicic et al., 2016). It produces a PCA plot computed from QC metrics, where cells corresponding to detected outliers are marked (Fig. 2b). Briefly, semi-robust principal components are computed from robustly sphered QC metrics data, using the pcout function from the mvoutlier R package (Filzmoser et al., 2008). These components are used to calculate distances for each cell, which are then used to compute weights for outlier detection. By default, the following QC metrics are used in this procedure: the percentage of counts from the top 100 features, the total number of features with detectable expression, the percentage of counts from control features, the number of detected feature controls, the log-scaled counts from endogenous features and log-transformed counts from feature controls. The user can specify QC metrics or other cell metadata variables to use for outlier detection with the ‘selected_variables’ argument to plotPCA. Detected outliers correspond to low-quality cells with abnormal library characteristics (e.g. low total counts and few expressed genes) that should be removed prior to downstream analysis. This automated approach is powerful but also somewhat opaque with respect to how outliers are defined, and so complements simpler filtering approaches that apply thresholds to particular QC metrics.

The plotQC function generates many types of plots useful for quality control, such as a plot to visualize the frequency of expression of features against their average expression level (Fig. 2c). Such plots are useful because scRNA-seq data are characterized by a high frequency of ‘dropout’ events, i.e. no observed expression (such as no read counts) in a particular cell for a gene that is actually expressed in that cell. Indeed, most genes will not have detectable expression in every cell. With plotQC, control features that should be present in each cell can be highlighted easily in the plot, allowing technical dropouts to be distinguished from biological heterogeneity of expression. Typical scRNA-seq datasets will show a broadly sigmoidal relationship between average expression level and frequency of expression across cells. This is consistent with expected behaviour where genes with greater average expression are more readily captured during library preparation and are detected at a greater frequency (Brennecke et al., 2013; Kim et al., 2015; Vallejos et al., 2015).

With plotQC, we can also produce a plot to visualize the most highly expressed features in the dataset (Fig. 2d). This provides a feature-centric overview of the dataset that visualizes the features with highest total expression across all cells, while also displaying the distribution of cell-level expression values for these features. It is common to see ERCC spike-ins (if used), mitochondrial and ribosomal genes among the highest expressed genes, while datasets consisting of healthy cells will also show high levels of constitutively expressed genes like ACTB. This plot allows the analyst to quickly check that the gene- or transcript-level quantification is behaving as expected, and to flag datasets where it is not.

Another important step in quality control is to identify variables (experimental factors or computed QC metrics) that drive variation in expression data across cells. The plotQC function provides a novel approach to identifying variables that have substantial explanatory power for many genes. For each variable in the phenoData slot of the SCESet object, we fit a linear model for each feature with only that variable as the explanatory variable. We then plot the distribution of the marginal R² values across all features for the variables with the most explanatory power for the dataset (Fig. 2e). The variables are ranked by median R² across features in the plot, allowing users to identify variables that may need to be considered during normalization or statistical modelling. The plotQC function can also assess the influence of variables of interest by plotting principal components of the expression matrix most strongly correlated with a variable of interest against that variable. For example, in the Case Study data, the first principal component is correlated with the C1 machine used to process the cell (Supplementary Fig. S2).

We also introduce the plotPhenoData function for convenient plotting of cell phenotype information (including QC metrics), and the plotFeatureData function for plotting feature information (see examples in the Supplementary Case Study). These methods will work not only on the SCESet class defined in scater, but also on any ExpressionSet object, providing sophisticated plotting functionality for many other Bioconductor packages and contexts.

The scater graphical user interface (GUI) provides convenient access to scater’s QC and visualization methods (Supplementary Figs S4–S6). This opens an interactive interface in a web browser that facilitates exploration of the data through QC plots and other intuitive visualizations. The GUI allows users of any background to easily examine the effects of changing multiple parameters, which can be helpful for quickly conducting exploratory data analysis. Useful settings can then be stored in R scripts to ensure that data analyses are reproducible.

In summary, scater provides a variety of novel and convenient methods to visualize an scRNA-seq dataset for QC. Low-quality cells and uninteresting genes can then be easily removed by filtering and subsetting the SCESet data structure prior to further analysis.

3.5 Data visualization

Dimensionality reduction techniques are necessary to convert high-dimensional expression data into low-dimensional representations for intuitive visualization of the relationships, similarities and differences between cells. To this end, scater provides convenient functions to apply a variety of dimensionality reduction procedures to the cells in an SCESet object. Functions include plotPCA, to perform a principal components analysis; plotTSNE, to perform t-distributed stochastic neighbour embedding (Van der Maaten and Hinton, 2008), which has been widely used for scRNA-seq data (Amir et al., 2013; Bendall et al., 2014; Macosko et al., 2015); plotDiffusionMap, to generate a diffusion map (Haghverdi et al., 2015) for visualizing differentiation processes; and plotMDS, to generate multi-dimensional scaling plots (Fig. 3a–c). The plotReducedDim function can also be used to plot any reduced-dimension representation of cells (e.g. an independent component analysis produced by monocle (Trapnell et al., 2013) or similar) that is stored in an SCESet object.

Open in a separate window

Fig. 3.

Reduced dimension representations of cells and gene expression plots with scater. Plots are shown using all genes (a–c) and cell cycle genes only (d–f) using PCA (a,d), t-SNE (b,e) and diffusion maps (c,f), where each point represents a cell. In the top row (a–c), points are coloured by patient of origin, sized by total features (number of genes with detectable expression) and the shape indicates the C1 machine used to process the cells. In the second row (d–f), points are coloured by the expression of CCND2 (a gene associated with the G1/S phase transition of the cell cycle) in each cell. Furthermore, with the plotExpression function, gene expression can be plotted against any cell metadata variables or the expression of another gene—here, expression for the CD86, IGH44 and IGHV4-34 genes in each cell is plotted against the patient of origin (g). The function automatically detects whether the x-axis variable is categorical or continuous and plots the data accordingly, with x-axis values ‘jittered’ to avoid excessive overplotting of points with the same x coordinate

By default, the PCA and t-SNE plots are produced using the features with the most variable expression across all cells. We focus on the most variable genes to highlight any heterogeneity in the data that might be driving interesting differences between cells. Alternatively, we can apply a priori knowledge to define a set of genes that are associated with a biological process of interest, and construct plots using only these features. For example, Scialdone et al. (2015) found that using prior knowledge to define feature sets is vital for exploring processes like the cell cycle, which can have substantial effects on single-cell expression measurements (Buettner et al., 2015). The subsetting and filtering methods for SCESet objects facilitate the generation of reduced-dimension plots for particular gene sets, in order to investigate certain effects in the data such as those due to the cell cycle (Fig. 3d–f).

The various types of reduced-dimension plots can be used to examine the structure of the cell population, including the formation of distinct subpopulations or the presence of continuous trajectories. Cell-level variables stored in the SCESet object can be used to define the shape, colour and size of points plotted, allowing more information to be conveniently incorporated into each plot (e.g. cells are coloured by CCND2 expression in Fig. 3d–f). The plotExpression function is also provided for plotting expression levels of a particular gene against any of the cell phenotype variables or the expression level of another feature (Fig. 3g). This allows the user to inspect the expression levels of a feature or set of features in full detail, rather than relying only on summary information and reduced-dimension plots where information is necessarily lost.

3.6 Data normalization and batch correction

Scaling normalization is typically required in RNA-seq data analysis to remove biases caused by differences in sequencing depth, capture efficiency or composition effects between samples. Frequently used methods for scaling normalization include the trimmed mean of M-values (Robinson and Oshlack, 2010), relative log-expression (Anders and Huber, 2010) and upper-quartile methods (Bullard et al., 2010), all of which are available for use in scater. In addition, scater is tightly integrated with the scran package, which implements a method utilizing cell pooling and deconvolution to compute size factors better suited to scRNA-seq data (Lun et al., 2016b). Lun et al. (2016a) also offers further discussion of the respective benefits and drawbacks of spike-in normalization and non-DE normalization.

After scaling normalization, further correction is typically required to ameliorate or remove batch effects. For example, in the case study dataset, cells from two patients were each processed on two C1 machines. Although C1 machine is not one of the most important explanatory variables on a per-gene level (Fig. 2e), this factor is correlated with the first principal component of the log-expression data (Fig. 2f). This effect cannot be removed by scaling normalization methods, which target cell-specific biases and are not sufficient for removing large-scale batch effects that vary on a gene-by-gene basis (Fig. 4a). Here we present two possibilities, both easily implemented in a scater workflow.

Open in a separate window

Fig. 4.

Normalization and batch correction with scater. Principal component analysis plots showing cell structure in the first two PCA dimensions using various normalization methods that can be easily applied in scater, including endogenous size-factor normalization using methods from the scran package (a); expression residuals after applying size-factor normalization and regressing out known, unwanted sources of variation (b); and removal of one hidden factor identified using the RUVs method from the RUV package (c). In all plots, the colour of points is determined by the patient from which cells were obtained, shape is determined by the C1 machine used to process the cells and size reflects the total number of genes with detectable expression in the cell

The C1 machine effect is known from the design of the experiment, so we can easily regress out this effect in scater. With the normaliseExprs function the user can supply a design matrix of variables to regress out of the expression values, and residuals from the linear model fit can be used as expression values for downstream analyses. For the dataset here, we fit a linear model to the scran normalized log-expression values with the C1 machine as an explanatory factor. (We also use the log-total counts from endogenous genes, percentage of counts from the top 100 most highly expressed genes and percentage of counts from control genes as additional covariates to control for these other unwanted technical effects.) We then use the residuals from the fitted model for further analyses (see Case Study in Supplementary Material). This approach successfully removes the C1 machine effect as a major source of variation between cells; the first principal component now separates the cells from the two patients, as expected (Fig. 4b). This approach needs to be used carefully as single-cell data often deviate from normal distributions, but in many cases, as here, it can successfully ameliorate large-scale known batch effects.

In addition to removing known batch effects, it can be important for large datasets to identify (potentially unknown) sources of unwanted variation (Grün and van Oudenaarden, 2015; Hicks et al., 2015; Leek et al., 2010). scater is compatible with existing methods such as svaseq (Leek and Storey, 2007; Leek, 2014) and RUVSeq (Risso et al., 2014) to identify and remove these unwanted sources of variation. Here, just removing the first latent variable identified by the RUVs method from RUVSeq is sufficient to remove the machine effect, as the PCA plot now separates cells by patient rather than C1 machine (Fig. 4c). More targeted applications of these methods can be used to remove specific effects, for example, by identifying latent factors from cell cycle genes to remove the cell cycle effect.

We emphasize that it is generally preferable to incorporate batch effects or latent variables into statistical models used for inference. Where this is not possible (e.g. for visualization), directly regressing out these uninteresting factors is required to obtain ‘corrected’ expression values for further analysis. Furthermore, a general risk of removing latent factors is that interesting biological variation may be removed along with the presumed unwanted variation. Users should therefore apply such methods with appropriate caution, particularly when an analysis aims to discover biological conditions, such as new cell types.

We thank Marco Salvetti for provision of the two samples and Zam Cader for processing the samples. We would also like to acknowledge Peter Donnelly and Oliver Stegle for support and helpful discussions.