Analysis of Large Single-Cell RNA-Seq Datasets in R/Bioconductor
Introduction
This workshop illustrates the latest Bioconductor infrastructure to analyze large single-cell datasets that may not entirely fit in memory.
We focus on the most common application in exploratory single-cell analysis, namely to find subpopulations of cells.
The proposed workflow consists of the following steps:
- Normalization;
- Dimensionality reduction;
- Clustering.
We will exploit a number of Bioconductor packages able to interact with the HDF5 on-disk data representation, freeing us of the need to load the full dataset in memory.
Note that this workshop is a high-level introduction to the analysis of large datasets from a user perspective. The focus is on packages that the user may directly use to analyze such data. The focus is NOT the direct interaction with infrastructural packages (such as the DelayedArray package) and is not aimed at developers wishing to build upon this infrastructure. If you are interested in these topics, see the workshop “Effectively using the DelayedArray framework to support the analysis of large data sets.”
Interacting with HDF5 files in R/Bioconductor
Here, we provide a brief introduction on the HDF5 infrastrucure in R/Bioconductor. See “Effectively using the DelayedArray framework to support the analysis of large data sets” for a more technical treatment of the topic.
At the low-level, the main interface between HDF5 and Bioconductor is implemented in the packages rhdf5, which provides read/write functionalities, Rhdf5lib, which provides C and C++ HDF5 libraries, and beachmat (Lun, Pagès, and Smith 2018), which provides a consistent C++ class interface for a variety of commonly used matrix types, including sparse and HDF5-backed matrices.
These packages are useful for developers that want to develop methods able to interact with HDF5 datasets. However, for most Bioconductor users interested in the analysis of single-cell data, the entry point is represented by the high-level class SingleCellExperiment (implemented in the SingleCellExperiment package ) and the lower level classes HDF5Matrix and DelayedMatrix, which can be stored in the assay slot of a SingleCellExperiment object. These matrix classes are implemented in the HDF5Array and DelayedArray packages, respectively. Once the data are stored in a SingleCellExperiment object with HDF5Matrix or DelayedMatrix as its assay, the packages scater (McCarthy et al. 2017), scran (Lun, Bach, and Marioni 2016), BiocSingular and mbkmeans can be seamlessly used.
The package DelayedMatrixStats deserves a special mention: it implements the rich API of the CRAN package matrixStats for HDF5Matrix and DelayedMatrix objects.
We invite the reader to find more details on all the mentioned packages in their relative vignettes. In the remainder of this workshop, we will use these methods to find cell sub-populations in a real datasets.
The SingleCellExperiment class
Add description of the class, with figure.
The dataset
To illustrate the analysis of large datasets, we use the largest of the Human Cell Atlas preview datasets available in the HCAData Bioconductor package.
The Human Cell Atlas is an international collaborative effort to map every cell type present in the human body (Regev et al. 2017). We focus on a set on the Census of Immune Cells dataset, specifically on the bone marrow subset (see https://preview.data.humancellatlas.org for additional details).
The dataset consists of 378,000 immune cells from the bone marrow of 8 healthy donors, processed with the 10X Genomics Chromium v2 platform at the Broad Institute. For each donor, 8 libraries were prepared, and 6000 cells were sequenced for each library.
Preprocessing
The data available in the package has already been preprocessed and the package contains a matrix of raw counts (in HDF5 format) of the top 6000 cell barcodes ranked by total UMI per barcode.
The raw counts were generated by Cell Ranger with GRCh38, standard 10X reference.
Note that:
- For each channel approximately 7000 cells were loaded, and the authors expect about 4000 non-empty droplets.
- Filtering of low quality barcodes is recommended.
More information on the preprocessing is available on the HCA website.
Loading the data in R
We use the HCAData package, which uses ExperimentHub to retrieve the data. Note that this is a time consuming step, because it needs to download the remotely hosted HDF5 matrix. DR: CONSIDER HOSTING SOMEWHERE JUST THE SUBSET?
start_time <- Sys.time()
library(HCAData)
library(ExperimentHub)
library(SingleCellExperiment)
eh <- ExperimentHub()
query(eh, "HCAData")## ExperimentHub with 6 records
## # snapshotDate(): 2019-06-20
## # $dataprovider: Human Cell Atlas
## # $species: Homo sapiens
## # $rdataclass: character
## # additional mcols(): taxonomyid, genome, description,
## # coordinate_1_based, maintainer, rdatadateadded, preparerclass,
## # tags, rdatapath, sourceurl, sourcetype
## # retrieve records with, e.g., 'object[["EH2047"]]'
##
## title
## EH2047 | Human Cell Atlas - Census of Immune Cells, Bone marrow, 'de...
## EH2048 | Human Cell Atlas - Census of Immune Cells, Bone marrow, sam...
## EH2049 | Human Cell Atlas - Census of Immune Cells, Bone marrow, gen...
## EH2050 | Human Cell Atlas - Census of Immune Cells, Umbilical cord b...
## EH2051 | Human Cell Atlas - Census of Immune Cells, Umbilical cord b...
## EH2052 | Human Cell Atlas - Census of Immune Cells, Umbilical cord b...sce <- HCAData("ica_bone_marrow")
sce## class: SingleCellExperiment
## dim: 33694 378000
## metadata(0):
## assays(1): counts
## rownames(33694): ENSG00000243485 ENSG00000237613 ...
## ENSG00000277475 ENSG00000268674
## rowData names(2): ID Symbol
## colnames(378000): MantonBM1_HiSeq_1-AAACCTGAGCAGGTCA-1
## MantonBM1_HiSeq_1-AAACCTGCACACTGCG-1 ...
## MantonBM8_HiSeq_8-TTTGTCATCTGCCAGG-1
## MantonBM8_HiSeq_8-TTTGTCATCTTGAGAC-1
## colData names(1): Barcode
## reducedDimNames(0):
## spikeNames(0):The object sce of class SingleCellExperiment contains the data in the assay counts, as well as a set of row and column data, accessible with the rowData() and colData() functions, respectively.
The assay counts is a DelayedMatrix with HDF5 backend that contains the gene expression data. It can be retrieve with either the counts() or assay() functions.
counts(sce)## <33694 x 378000> DelayedMatrix object of type "integer":
## MantonBM1_HiSeq_1-AAACCTGAGCAGGTCA-1 ...
## ENSG00000243485 0 .
## ENSG00000237613 0 .
## ENSG00000186092 0 .
## ENSG00000238009 0 .
## ENSG00000239945 0 .
## ... . .
## ENSG00000277856 0 .
## ENSG00000275063 0 .
## ENSG00000271254 0 .
## ENSG00000277475 0 .
## ENSG00000268674 0 .
## MantonBM8_HiSeq_8-TTTGTCATCTTGAGAC-1
## ENSG00000243485 0
## ENSG00000237613 0
## ENSG00000186092 0
## ENSG00000238009 0
## ENSG00000239945 0
## ... .
## ENSG00000277856 0
## ENSG00000275063 0
## ENSG00000271254 0
## ENSG00000277475 0
## ENSG00000268674 0To check how where data is stored, we can use the seed() function, which gives valuable information on the geometry of the HDF5 file.
seed(counts(sce))## An object of class "HDF5ArraySeed"
## Slot "filepath":
## [1] "/home/rstudio/.cache/ExperimentHub/325c13401453_2047"
##
## Slot "name":
## [1] "counts"
##
## Slot "dim":
## [1] 33694 378000
##
## Slot "first_val":
## [1] 0
##
## Slot "chunkdim":
## [1] 184 2060Library information
Unfortunately, the column data does not include any metadata on the libraries, although some is available at the HCA website.
We can extract information on the donors and libraries by parsing the barcodes.
library(stringr)
cell_info <- str_split(colData(sce)$Barcode, "_", simplify = TRUE)
donor <- cell_info[,1]
lib <- paste0(donor, "_", str_split(cell_info[,3], "-", simplify = TRUE)[,1])
colData(sce)$donor <- donor
colData(sce)$lib <- libWe add donor and library information to the column data. This will be useful for subsetting.
Subsetting
To work with a small enough dataset, we will randomly subsample the data to have only 50 cells per library. This will allow us to run the workshop in a reasonable time, while preserving the batch structure typical of a large-scale experiment.
Note that it is sufficient to add eval=FALSE to the next code chunk to run the workflow on the full dataset.
library(dplyr)
set.seed(50)
colData(sce) %>%
as.data.frame() %>%
group_by(lib) %>%
sample_n(50) %>%
pull(Barcode) -> cell_idx
sce <- sce[, cell_idx]
sce## class: SingleCellExperiment
## dim: 33694 3150
## metadata(0):
## assays(1): counts
## rownames(33694): ENSG00000243485 ENSG00000237613 ...
## ENSG00000277475 ENSG00000268674
## rowData names(2): ID Symbol
## colnames(3150): MantonBM1_HiSeq_1-TGCCCTATCGGCGCTA-1
## MantonBM1_HiSeq_1-GGGAATGAGTAACCCT-1 ...
## MantonBM8_HiSeq_8-GGACATTTCCTACAGA-1
## MantonBM8_HiSeq_8-TCAATCTTCGTTGCCT-1
## colData names(3): Barcode donor lib
## reducedDimNames(0):
## spikeNames(0):counts(sce)## <33694 x 3150> DelayedMatrix object of type "integer":
## MantonBM1_HiSeq_1-TGCCCTATCGGCGCTA-1 ...
## ENSG00000243485 0 .
## ENSG00000237613 0 .
## ENSG00000186092 0 .
## ENSG00000238009 0 .
## ENSG00000239945 0 .
## ... . .
## ENSG00000277856 0 .
## ENSG00000275063 0 .
## ENSG00000271254 0 .
## ENSG00000277475 0 .
## ENSG00000268674 0 .
## MantonBM8_HiSeq_8-TCAATCTTCGTTGCCT-1
## ENSG00000243485 0
## ENSG00000237613 0
## ENSG00000186092 0
## ENSG00000238009 0
## ENSG00000239945 0
## ... .
## ENSG00000277856 0
## ENSG00000275063 0
## ENSG00000271254 0
## ENSG00000277475 0
## ENSG00000268674 0Note that because we have not actually made any computations, the original HDF5 file is untouched.
seed(counts(sce))## An object of class "HDF5ArraySeed"
## Slot "filepath":
## [1] "/home/rstudio/.cache/ExperimentHub/325c13401453_2047"
##
## Slot "name":
## [1] "counts"
##
## Slot "dim":
## [1] 33694 378000
##
## Slot "first_val":
## [1] 0
##
## Slot "chunkdim":
## [1] 184 2060To avoid reading from this huge file, it is best to save the subset on a new, smaller HDF5 file. This can be done with the realize() function.
counts(sce) <- realize(counts(sce), BACKEND = "HDF5Array")
seed(counts(sce))## An object of class "HDF5ArraySeed"
## Slot "filepath":
## [1] "/tmp/RtmpkXNUiN/HDF5Array_dump/auto00001.h5"
##
## Slot "name":
## [1] "/HDF5ArrayAUTO00001"
##
## Slot "dim":
## [1] 33694 3150
##
## Slot "first_val":
## [1] 0
##
## Slot "chunkdim":
## [1] 3270 305This is a time consuming step, but will make the subsequent calls much faster.
Parallel computing
For the workshop, we run the workflow in serial mode. When running the workflow on large datasets, we recommend running the workflow in parallel, by using the BiocParallel package. See the BiocParallel vignette for more information on parallel computing.
library(BiocParallel)
register(SerialParam())
register(MulticoreParam(6)) ## remove this!Quality control
Filtering of low-quality samples
First, we use the scater package to compute a set of QC measures and filter out the low-quality samples.
library(scater)
sce <- calculateQCMetrics(sce,
feature_controls=list(Mito=grep("^MT", rowData(sce)$Symbol)))We remove cells with low UMI counts and high proportion of mitocondrial reads, using it as a proxy for cell damage.
high_mito <- isOutlier(sce$pct_counts_Mito, nmads=3, type="higher")
low_counts <- sce$total_counts < 100
table((!high_mito) & (!low_counts))##
## FALSE TRUE
## 412 2738sce <- sce[, (!high_mito) & (!low_counts)]
sce## class: SingleCellExperiment
## dim: 33694 2738
## metadata(0):
## assays(1): counts
## rownames(33694): ENSG00000243485 ENSG00000237613 ...
## ENSG00000277475 ENSG00000268674
## rowData names(10): ID Symbol ... total_counts log10_total_counts
## colnames(2738): MantonBM1_HiSeq_1-TGCCCTATCGGCGCTA-1
## MantonBM1_HiSeq_1-GGGAATGAGTAACCCT-1 ...
## MantonBM8_HiSeq_8-CTAGCCTGTTTACTCT-1
## MantonBM8_HiSeq_8-GGACATTTCCTACAGA-1
## colData names(39): Barcode donor ...
## pct_counts_in_top_200_features_Mito
## pct_counts_in_top_500_features_Mito
## reducedDimNames(0):
## spikeNames(0):Removal of lowly expressed genes
In addition, we keep only those genes that have at least 1 UMI in at least 5% of the data. These threshold are dataset-specific and may need to be taylored to specific applications.
num_reads <- 1
num_cells <- 0.05*ncol(sce)
keep <- which(DelayedArray::rowSums(counts(sce) >= num_reads ) >= num_cells)
sce <- sce[keep,]
sce## class: SingleCellExperiment
## dim: 4442 2738
## metadata(0):
## assays(1): counts
## rownames(4442): ENSG00000279457 ENSG00000188976 ...
## ENSG00000198695 ENSG00000198727
## rowData names(10): ID Symbol ... total_counts log10_total_counts
## colnames(2738): MantonBM1_HiSeq_1-TGCCCTATCGGCGCTA-1
## MantonBM1_HiSeq_1-GGGAATGAGTAACCCT-1 ...
## MantonBM8_HiSeq_8-CTAGCCTGTTTACTCT-1
## MantonBM8_HiSeq_8-GGACATTTCCTACAGA-1
## colData names(39): Barcode donor ...
## pct_counts_in_top_200_features_Mito
## pct_counts_in_top_500_features_Mito
## reducedDimNames(0):
## spikeNames(0):Again, it is worth realize()-ing the matrix to speed up calculations.
counts(sce) <- realize(counts(sce), BACKEND = "HDF5Array")
seed(counts(sce))## An object of class "HDF5ArraySeed"
## Slot "filepath":
## [1] "/tmp/RtmpkXNUiN/HDF5Array_dump/auto00002.h5"
##
## Slot "name":
## [1] "/HDF5ArrayAUTO00002"
##
## Slot "dim":
## [1] 4442 2738
##
## Slot "first_val":
## [1] 0
##
## Slot "chunkdim":
## [1] 1273 785Normalization
Normalization is a crucial step in the preprocessing of the results. Here, we use the scran package to compute size factors that we will use to compute the normalized log-expression values.
It has been shown that the scran method works best if the size factors are computed within roughly homogeneous cell populations; hence, it is beneficial to run a quick clustering on the raw data to compute better size factors. This ensures that we do not pool cells that are very different. Note that this is not the final clustering to identify cell sub-populations.
In order to carry out this clustering, we use a mini-batch k-means algorithm, implemented in the mbkmeans package. We will talk more about the details of this method in the Clustering section.
DR: CONSIDER NORMALIZING EACH BATCH INDEPENDENTLY AND THEN MERGING THEM WITH MNN.
library(scran)
library(mbkmeans)
set.seed(1000)
clusters <- mbkmeans(counts(sce), clusters=10, batch_size = 100)
table(clusters$Clusters)##
## 1 2 3 4 5 6 7 8 9 10
## 3 8 12 430 667 6 799 715 64 34sce <- computeSumFactors(sce, min.mean=0.1, cluster=clusters$Clusters)It is useful to check whether the size factors are correlated with the total number of reads per cell.
plot(sce$total_counts, sizeFactors(sce), log="xy", xlab="Total reads", ylab="scran size factors")
Finally, we compute normalized log-expression values with the normalize() function from the scater package.
sce <- scater::normalize(sce)Note that the log-normalized data are stored in the logcounts assay of the object. Since the counts assay is a DelayedMatrix and we have only one set of size factors in the object, also the normalized data are stored as a DelayedMatrix.
logcounts(sce)## <4442 x 2738> DelayedMatrix object of type "double":
## MantonBM1_HiSeq_1-TGCCCTATCGGCGCTA-1 ...
## ENSG00000279457 0 .
## ENSG00000188976 0 .
## ENSG00000187608 0 .
## ENSG00000078808 0 .
## ENSG00000160087 0 .
## ... . .
## ENSG00000212907 0.000000 .
## ENSG00000198886 3.996768 .
## ENSG00000198786 2.687402 .
## ENSG00000198695 0.000000 .
## ENSG00000198727 3.395429 .
## MantonBM8_HiSeq_8-GGACATTTCCTACAGA-1
## ENSG00000279457 0.0000000
## ENSG00000188976 0.2639181
## ENSG00000187608 0.0000000
## ENSG00000078808 0.0000000
## ENSG00000160087 0.0000000
## ... .
## ENSG00000212907 0.4869413
## ENSG00000198886 4.0753686
## ENSG00000198786 1.7692736
## ENSG00000198695 0.0000000
## ENSG00000198727 3.7273595This allows us to store in memory only the colData and rowData, resulting in a fairly small object.
library(pryr)
object_size(sce)## 2.71 MBDimensionality reduction
library(BiocSingular)
library(DelayedMatrixStats)
## find most variable genes
vars <- DelayedMatrixStats::rowVars(logcounts(sce))
names(vars) <- rownames(sce)
vars <- sort(vars, decreasing = TRUE)
for_pca <- t(logcounts(sce)[names(vars)[1:1000],])
pca <- BiocSingular::runPCA(for_pca, rank = 30,
scale = TRUE,
BSPARAM = RandomParam(deferred = FALSE))
reducedDim(sce, "PCA") <- pca$x
sce## class: SingleCellExperiment
## dim: 4442 2738
## metadata(1): log.exprs.offset
## assays(2): counts logcounts
## rownames(4442): ENSG00000279457 ENSG00000188976 ...
## ENSG00000198695 ENSG00000198727
## rowData names(10): ID Symbol ... total_counts log10_total_counts
## colnames(2738): MantonBM1_HiSeq_1-TGCCCTATCGGCGCTA-1
## MantonBM1_HiSeq_1-GGGAATGAGTAACCCT-1 ...
## MantonBM8_HiSeq_8-CTAGCCTGTTTACTCT-1
## MantonBM8_HiSeq_8-GGACATTTCCTACAGA-1
## colData names(39): Barcode donor ...
## pct_counts_in_top_200_features_Mito
## pct_counts_in_top_500_features_Mito
## reducedDimNames(1): PCA
## spikeNames(0):plotPCA(sce, colour_by = "donor")
Clustering with Mini-batch k-means
Choose the number of clusters
wcss <- lapply(10:20, function(k) {
cl <- mbkmeans(sce, reduceMethod = "PCA", clusters = k,
batch_size = 100, num_init=10, max_iters=100,
calc_wcss = TRUE)
})plot(10:20, sapply(wcss, function(x) sum(x$WCSS_per_cluster)),
type="b",
xlab="Number of clusters K",
ylab="Total within-clusters sum of squares")
clusters2 <- mbkmeans(sce, reduceMethod = "PCA", clusters=16,
batch_size=500,
num_init = 10, max_iters = 100)
sce$cluster <- as.factor(clusters2$Clusters)Visualization
sce <- runTSNE(sce, use_dimred = "PCA",
external_neighbors=TRUE,
BNPARAM = BiocNeighbors::AnnoyParam(),
nthreads = 1)
plotTSNE(sce, colour_by = "cluster")
plotTSNE(sce, colour_by = "donor")
plotTSNE(sce, colour_by = "ENSG00000090382") + ggtitle("LYZ (CD14+)")
plotTSNE(sce, colour_by = "ENSG00000156738") + ggtitle("MS4A1 (B-cells)")
plotTSNE(sce, colour_by = "ENSG00000105374") + ggtitle("NKG7 (NK-cells)")
sce <- runUMAP(sce, use_dimred = "PCA",
external_neighbors=TRUE,
BNPARAM = BiocNeighbors::AnnoyParam())
plotUMAP(sce, colour_by = "cluster")
plotUMAP(sce, colour_by = "donor")
DR: CONSIDER ADDING BATCH EFFECT CORRECTION?
Time
print(Sys.time() - start_time)## Time difference of 5.228855 minsSession Info
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 9 (stretch)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [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] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] DelayedMatrixStats_1.7.1 BiocSingular_1.1.3
## [3] pryr_0.1.4 mbkmeans_1.1.0
## [5] scran_1.13.7 scater_1.13.7
## [7] ggplot2_3.2.0 HDF5Array_1.13.2
## [9] dplyr_0.8.1 stringr_1.4.0
## [11] rhdf5_2.29.0 ExperimentHub_1.11.1
## [13] AnnotationHub_2.17.3 BiocFileCache_1.9.1
## [15] dbplyr_1.4.2 HCAData_1.1.1
## [17] SingleCellExperiment_1.7.0 SummarizedExperiment_1.15.3
## [19] DelayedArray_0.11.2 BiocParallel_1.19.0
## [21] matrixStats_0.54.0 Biobase_2.45.0
## [23] GenomicRanges_1.37.12 GenomeInfoDb_1.21.1
## [25] IRanges_2.19.10 S4Vectors_0.23.13
## [27] BiocGenerics_0.31.4 BiocManager_1.30.4
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 ggbeeswarm_0.6.0
## [3] colorspace_1.4-1 dynamicTreeCut_1.63-1
## [5] benchmarkme_1.0.0 XVector_0.25.0
## [7] BiocNeighbors_1.3.2 bit64_0.9-7
## [9] interactiveDisplayBase_1.23.0 AnnotationDbi_1.47.0
## [11] ClusterR_1.1.9 codetools_0.2-16
## [13] splines_3.6.0 doParallel_1.0.14
## [15] knitr_1.23 ade4_1.7-13
## [17] cluster_2.0.8 uwot_0.1.3
## [19] FD_1.0-12 shiny_1.3.2
## [21] compiler_3.6.0 httr_1.4.0
## [23] dqrng_0.2.1 assertthat_0.2.1
## [25] Matrix_1.2-17 lazyeval_0.2.2
## [27] limma_3.41.5 later_0.8.0
## [29] htmltools_0.3.6 tools_3.6.0
## [31] gmp_0.5-13.5 rsvd_1.0.1
## [33] igraph_1.2.4.1 gtable_0.3.0
## [35] glue_1.3.1 GenomeInfoDbData_1.2.1
## [37] rappdirs_0.3.1 Rcpp_1.0.1
## [39] ape_5.3 nlme_3.1-139
## [41] blogdown_0.13.1 iterators_1.0.10
## [43] xfun_0.7 beachmat_2.1.1
## [45] mime_0.7 irlba_2.3.3
## [47] gtools_3.8.1 statmod_1.4.32
## [49] edgeR_3.27.5 zlibbioc_1.31.0
## [51] MASS_7.3-51.4 scales_1.0.0
## [53] promises_1.0.1 yaml_2.2.0
## [55] curl_3.3 memoise_1.1.0
## [57] gridExtra_2.3 stringi_1.4.3
## [59] RSQLite_2.1.1 foreach_1.4.4
## [61] permute_0.9-5 benchmarkmeData_1.0.1
## [63] geometry_0.4.1 rlang_0.3.4
## [65] pkgconfig_2.0.2 bitops_1.0-6
## [67] evaluate_0.14 lattice_0.20-38
## [69] purrr_0.3.2 Rhdf5lib_1.7.2
## [71] labeling_0.3 cowplot_0.9.4
## [73] bit_1.1-14 tidyselect_0.2.5
## [75] magrittr_1.5 bookdown_0.11.1
## [77] R6_2.4.0 DBI_1.0.0
## [79] pillar_1.4.1 withr_2.1.2
## [81] mgcv_1.8-28 abind_1.4-5
## [83] RCurl_1.95-4.12 tibble_2.1.3
## [85] crayon_1.3.4 rmarkdown_1.13
## [87] viridis_0.5.1 locfit_1.5-9.1
## [89] grid_3.6.0 blob_1.1.1
## [91] vegan_2.5-5 digest_0.6.19
## [93] xtable_1.8-4 httpuv_1.5.1
## [95] RcppParallel_4.4.3 munsell_0.5.0
## [97] beeswarm_0.2.3 viridisLite_0.3.0
## [99] vipor_0.4.5 magic_1.5-9References
Lun, Aaron TL, Karsten Bach, and John C Marioni. 2016. “Pooling Across Cells to Normalize Single-Cell Rna Sequencing Data with Many Zero Counts.” Genome Biology 17 (1). BioMed Central: 75.
Lun, Aaron TL, Hervé Pagès, and Mike L Smith. 2018. “Beachmat: A Bioconductor C++ Api for Accessing High-Throughput Biological Data from a Variety of R Matrix Types.” PLoS Computational Biology 14 (5). Public Library of Science: e1006135.
McCarthy, Davis J, Kieran R Campbell, Aaron TL Lun, and Quin F Wills. 2017. “Scater: Pre-Processing, Quality Control, Normalization and Visualization of Single-Cell Rna-Seq Data in R.” Bioinformatics 33 (8). Oxford University Press: 1179–86.
Regev, Aviv, Sarah A Teichmann, Eric S Lander, Ido Amit, Christophe Benoist, Ewan Birney, Bernd Bodenmiller, et al. 2017. “Science Forum: The Human Cell Atlas.” eLife 6. eLife Sciences Publications Limited: e27041.