scater: Single-cell analysis toolkit for expression in Rscater 1.8.4
This document gives an introduction to and overview of the quality control functionality of the scater package. scater contains tools to help with the analysis of single-cell transcriptomic data, focusing on RNA-seq data. The package features:
SingleCellExperiment class as a data container for interoperability with a wide range of other Bioconductor packages;SingleCellExperiment objectWe assume that you have a matrix containing expression count data summarised at the level of some features (gene, exon, region, etc.).
First, we create a SingleCellExperiment object containing the data, as demonstrated below with some example data ("sc_example_counts") and metadata ("sc_example_cell_info"):
Rows of the object correspond to features, while columns correspond to samples, i.e., cells in the context of single-cell ’omics data.
library(scater)
data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
assays = list(counts = sc_example_counts),
colData = sc_example_cell_info
)
example_sce## class: SingleCellExperiment
## dim: 2000 40
## metadata(0):
## assays(1): counts
## rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
## rowData names(0):
## colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
## colData names(4): Cell Mutation_Status Cell_Cycle Treatment
## reducedDimNames(0):
## spikeNames(0):We usually expect (raw) count data to be labelled as "counts" in the assays, which can be easily retrieved with the counts accessor.
Getters and setters are also provided for exprs, tpm, cpm, fpkm and versions of these with the prefix norm_.
str(counts(example_sce))Row and column-level metadata are easily accessed (or modified) as shown below.
There are also dedicated getters and setters for spike-in specifiers (isSpike); size factor values (sizeFactors); and reduced dimensionality results (reducedDim).
example_sce$whee <- sample(LETTERS, ncol(example_sce), replace=TRUE)
colData(example_sce)## DataFrame with 40 rows and 5 columns
## Cell Mutation_Status Cell_Cycle Treatment whee
## <character> <character> <character> <character> <character>
## Cell_001 Cell_001 positive S treat1 U
## Cell_002 Cell_002 positive G0 treat1 W
## Cell_003 Cell_003 negative G1 treat1 J
## Cell_004 Cell_004 negative S treat1 E
## Cell_005 Cell_005 negative G1 treat2 Z
## ... ... ... ... ... ...
## Cell_036 Cell_036 negative G0 treat1 P
## Cell_037 Cell_037 negative G0 treat1 I
## Cell_038 Cell_038 negative G0 treat2 F
## Cell_039 Cell_039 negative G1 treat1 P
## Cell_040 Cell_040 negative G0 treat2 ProwData(example_sce)$stuff <- runif(nrow(example_sce))
rowData(example_sce)## DataFrame with 2000 rows and 1 column
## stuff
## <numeric>
## 1 0.0673663057386875
## 2 0.222470983397216
## 3 0.947641894221306
## 4 0.374451905954629
## 5 0.593446868704632
## ... ...
## 1996 0.474242345895618
## 1997 0.192777156364173
## 1998 0.151339191012084
## 1999 0.400677960831672
## 2000 0.937458637170494Subsetting is very convenient with this class, as both data and metadata are processed in a synchronized manner. For example, we can filter out features (genes) that are not expressed in any cells:
keep_feature <- rowSums(counts(example_sce) > 0) > 0
example_sce <- example_sce[keep_feature,]More details about the SingleCellExperiment class can be found in the documentation for SingleCellExperiment package.
We calculate counts-per-million using the aptly-named calculateCPM function.
The output is most appropriately stored as an assay named "cpm" in the assays of the SingleCellExperiment object.
cpm(example_sce) <- calculateCPM(example_sce)Another option is to use the normalize function, which calculates log2-transformed normalized expression values.
This is done by dividing each count by its size factor (or scaled library size, if no size factors are defined), adding a pseudo-count and log-transforming.
The resulting values can be interpreted on the same scale as log-transformed counts, and are stored in "logcounts".
example_sce <- normalize(example_sce)
assayNames(example_sce)## [1] "counts" "cpm" "logcounts"Note that exprs is a synonym for logcounts when accessing or setting data.
This is done for backwards compatibility with older verions of scater.
identical(exprs(example_sce), logcounts(example_sce))## [1] TRUEOf course, users can construct any arbitrary matrix of the same dimensions as the count matrix and store it as an assay.
assay(example_sce, "is_expr") <- counts(example_sce)>0Count matrices stored as CSV files or equivalent can be easily read into R session using read.table from utils or fread from the data.table package.
It is advisable to coerce the resulting object into a matrix before storing it in a SingleCellExperiment object.
Data from 10X Genomics experiments can be read in using the read10xCounts function from the DropletUtils package.
This will automatically generate a SingleCellExperiment with a sparse matrix, see the documentation for more details.
scater also provides wrapper functions to kallisto and Salmon for easy and extremely fast quantification of transcript abundance from within R.
An example of how to use these functions is shown below:
# Example if in the kallisto/test directory
setwd("/home/davis/kallisto/test")
kallisto_log <- runKallisto("targets.txt", "transcripts.idx", single_end=FALSE,
output_prefix="output", verbose=TRUE, n_bootstrap_samples=10)
sce_test <- readKallistoResults(kallisto_log, read_h5=TRUE)
sce_testSCESet classAs of July 2017, scater has switched from the SCESet class previously defined within the package to the more widely applicable SingleCellExperiment class.
From Bioconductor 3.6 (October 2017), the release version of scater will use SingleCellExperiment.
SingleCellExperiment is a more modern and robust class that provides a common data structure used by many single-cell Bioconductor packages.
Advantages include support for sparse data matrices and the capability for on-disk storage of data to minimise memory usage for large single-cell datasets.
It should be straight-forward to convert existing scripts based on SCESet objects to SingleCellExperiment objects, with key changes outlined immediately below.
toSingleCellExperiment and updateSCESet (for backwards compatibility) can be used to convert an old SCESet object to a SingleCellExperiment object;SingleCellExperiment object with the function SingleCellExperiment (actually less fiddly than creating a new SCESet);scater functions have been refactored to take SingleCellExperiment objects, so once data is in a SingleCellExperiment object, the user experience is almost identical to that with the SCESet class.Users may need to be aware of the following when updating their own scripts:
colnames function (instead of sampleNames or cellNames for an SCESet object);rownames function (instead of featureNames);phenoData in an SCESet, corresponds to colData in a SingleCellExperiment object and is accessed/assigned with the colData function (this replaces the pData function);$ operator (e.g. sce$total_counts);featureData in an SCESet, corresponds to rowData in a SingleCellExperiment object and is accessed/assigned with the rowData function (this replaces the fData function);plotScater, which produces a cumulative expression, overview plot, replaces
the generic plot function for SCESet objects.sessionInfo()## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.7-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] scater_1.8.4 SingleCellExperiment_1.2.0
## [3] SummarizedExperiment_1.10.1 DelayedArray_0.6.4
## [5] BiocParallel_1.14.2 matrixStats_0.54.0
## [7] GenomicRanges_1.32.6 GenomeInfoDb_1.16.0
## [9] IRanges_2.14.10 S4Vectors_0.18.3
## [11] Biobase_2.40.0 BiocGenerics_0.26.0
## [13] ggplot2_3.0.0 knitr_1.20
## [15] BiocStyle_2.8.2
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.13 ggbeeswarm_0.6.0
## [3] colorspace_1.3-2 RcppEigen_0.3.3.4.0
## [5] rjson_0.2.20 class_7.3-14
## [7] rio_0.5.10 rprojroot_1.3-2
## [9] XVector_0.20.0 proxy_0.4-22
## [11] tximport_1.8.0 robustbase_0.93-2
## [13] shinydashboard_0.7.0 shiny_1.1.0
## [15] compiler_3.5.1 backports_1.1.2
## [17] assertthat_0.2.0 Matrix_1.2-14
## [19] lazyeval_0.2.1 limma_3.36.2
## [21] later_0.7.3 htmltools_0.3.6
## [23] tools_3.5.1 bindrcpp_0.2.2
## [25] igraph_1.2.2 gtable_0.2.0
## [27] glue_1.3.0 GenomeInfoDbData_1.1.0
## [29] reshape2_1.4.3 dplyr_0.7.6
## [31] ggthemes_4.0.0 Rcpp_0.12.18
## [33] carData_3.0-1 cellranger_1.1.0
## [35] DelayedMatrixStats_1.2.0 lmtest_0.9-36
## [37] xfun_0.3 laeken_0.4.6
## [39] stringr_1.3.1 openxlsx_4.1.0
## [41] mime_0.5 edgeR_3.22.3
## [43] DEoptimR_1.0-8 zlibbioc_1.26.0
## [45] MASS_7.3-50 zoo_1.8-3
## [47] scales_1.0.0 VIM_4.7.0
## [49] hms_0.4.2 promises_1.0.1
## [51] rhdf5_2.24.0 yaml_2.2.0
## [53] curl_3.2 gridExtra_2.3
## [55] stringi_1.2.4 e1071_1.7-0
## [57] destiny_2.10.2 TTR_0.23-3
## [59] boot_1.3-20 zip_1.0.0
## [61] rlang_0.2.1 pkgconfig_2.0.1
## [63] bitops_1.0-6 evaluate_0.11
## [65] lattice_0.20-35 purrr_0.2.5
## [67] Rhdf5lib_1.2.1 bindr_0.1.1
## [69] labeling_0.3 cowplot_0.9.3
## [71] tidyselect_0.2.4 plyr_1.8.4
## [73] magrittr_1.5 bookdown_0.7
## [75] R6_2.2.2 pillar_1.3.0
## [77] haven_1.1.2 foreign_0.8-71
## [79] withr_2.1.2 xts_0.11-0
## [81] scatterplot3d_0.3-41 abind_1.4-5
## [83] RCurl_1.95-4.11 sp_1.3-1
## [85] nnet_7.3-12 tibble_1.4.2
## [87] crayon_1.3.4 car_3.0-0
## [89] rmarkdown_1.10 viridis_0.5.1
## [91] locfit_1.5-9.1 grid_3.5.1
## [93] readxl_1.1.0 data.table_1.11.4
## [95] forcats_0.3.0 vcd_1.4-4
## [97] digest_0.6.15 xtable_1.8-2
## [99] httpuv_1.4.5 munsell_0.5.0
## [101] beeswarm_0.2.3 viridisLite_0.3.0
## [103] smoother_1.1 vipor_0.4.5