title	author	date	package	abstract	vignette	output
Accessing TENxGenomics data in _R_ / _Bioconductor_	Martin Morgan	`r doc_date()`	`r pkg_ver('TENxGenomics')`	`r packageDescription('TENxGenomics')$Description`	%\VignetteIndexEntry{Accessing TENxGenomics data in R / Bioconductor} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8}	BiocStyle::html_document

knitr::opts_chunk$set(
    eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")),
    cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE"))
)
suppressPackageStartupMessages({
    library(TENxGenomics)
    library(BiocFileCache)
    library(SummarizedExperiment)
    library(Rtsne)
})

Setup

This vignette requires the TENxGenomics package, available from github.

biocLite("mtmorgan/TENxGenomics")
library(TENxGenomics)

The vignette uses large datasets made available from 10xGenomics. We store these in a convenient location using BiocFileCache.

library(BiocFileCache)
bfc <- BiocFileCache()

oneM <- paste0(
    "https://s3-us-west-2.amazonaws.com/10x.files/",
    "samples/cell/1M_neurons/",
    "1M_neurons_filtered_gene_bc_matrices_h5.h5"
)
path <- bfcrpath(bfc, oneM)

Discovery and subsetting

The 10x data are 'hdf5' format files. Discover basic information about the data set using the TENxGenomics() constructor.

tenx <- TENxGenomics(path)
tenx

The returned object is a light-weight 'view' into the file. The view has matrix-like semantics, with methods dim() (implicitly, nrow(), ncol()), dimnames() (rownames() and colnames()), and [. The latter is useful to easily subset the very large data to a more useful size. Subsetting supports numeric, character, and logical vectors.

tenx[, sample(ncol(tenx), 1000)]
colnames(tenx[, sample(ncol(tenx), 3)])

Input

A useful strategy when working with large data is to input portions of the data. This allows, for instance, management of overall memory use when exploiting multiple computational cores. On typical computers it might be reasonable to input on the order of 10k samples at a time.

Simple

Use as.matrix() (dense matrix) or as.dgCMatrix() (sparse matrix representation) to read a subset of the actual data in to R.

onek <- as.matrix(tenx[, 1:1000])
class(onek)
dim(onek)
onek[1:10, 1:5]

Input is quickest when the columns are sequential, but one can also input random rows and columns. This is reasonably quick for samples up to about 1k.

as.matrix(tenx[sample(nrow(tenx), 5), sample(ncol(tenx), 3)])

Using a TENxMatrix object

An alternative to creating TENxGenomics object tenx is to wrap the 10xGenomics data in a TENxMatrix object.

tenxmat <- TENxMatrix(path)

The TENxMatrix class extends the DelayedArray class defined in the DelayedArray package so all the operations available on DelayedArray objects work on TENxMatrix objects. See ?DelayedArray for more information.

Rich

It is often helpful to place raw count data such as that returned by as.matrix() or as.dgCMatrix() into experimental context, e.g., the cell, library, and mouse from which the information has been derived. The SummarizedExperiment package and class is the standard Bioconductor container for this type of representation.

Here we create a SummarizedExperiment around the TENxGenomics representation. The object infers information (as described on ?tenxSummarizedExperiment) about the library and mouse brain used for each sample. We use this to identify 100 random cells from mouse "A", and 100 random cells from mouse "B".

tenxse <- tenxSummarizedExperiment(path)
colData(tenxse)
n <- 100
samples <- as.vector(vapply(
    split(tenxse$Barcode, tenxse$Mouse),
    sample, character(n), n
))

We then instantiate the data as a matrix in a SummarizedExperiment, either directly from the file path, or from a TENxGenomics instance.

library(SummarizedExperiment)
se <- matrixSummarizedExperiment(path, j = samples)
se
table(se$Mouse)

Iterative

Simple or rich input is useful when wishing to work with a portion of the data that fits in memory, especially during exploratory phases of analysis. Processing the whole file requires some kind of iterative approach because, like all programming lagauges, it makes little sense to read very large volumes of data into main memory. The tenxreduce() function visits the entire hdf5 file, return column-oriented slices filtered through the rows and columns present in the TENxGenomics argument.

Here we use a smaller data set for illustrative purposes

twentyK <- paste0(
    "https://s3-us-west-2.amazonaws.com/10x.files/",
    "samples/cell/1M_neurons/",
    "1M_neurons_neuron20k.h5"
)
path <- bfcrpath(bfc, twentyK)
tenx <- TENxGenomics(path)
tenx

The tenxiterate() function takes a TENxGenomics instance and a function FUN(). FUN() accepts at least one argument, e.g., x. FUN(x, ...) is called on successive chunks of the hdf5 file. The argument x is a list, with elements containing the row index (x$ridx), column index (x$cidx), and read count (x$value) of a slice of the hdf5 data. FUN() peforms arbitrary transformations on the data, and the result is accumulated across chunks. The function is implemented on top of BiocParallel::bpiterate(), so supports parallel processing. The following processes the data in chunks, calculating the total number of aligned reads.

BiocParallel::register(
    BiocParallel::MulticoreParam(progressbar = FALSE)
)
result <- tenxiterate(tenx, function(x) sum(x$value))  # reads per chunk
sum(unlist(result))                                    # reads total

The following summarizes the row and column margins, with n the number of non-zero cells and sum the number of reads per row or column. The chunks are 'sparse' representations, with continguous columns, so efficient processing takes different stratgies. Some care is also taken to reduce (though not minimize) the size of data returned by the function, for better performance when evaluated in a parallel context.

margin.summary <- function(x, nrow) {
    ## > str(x)
    ## List of 3
    ##  $ ridx : num [1:20548381] 8 9 17 39 52 63 118 119 123 182 ...
    ##  $ cidx : num [1:20548381] 1 1 1 1 1 1 1 1 1 1 ...
    ##  $ value: int [1:20548381] 1 1 2 2 1 7 2 1 1 1 ...

    ## rows: summarize all rows, whether in current sample or not.
    ridx <- structure(                  # quick 'factor'
        x$ridx, .Label=as.character(seq_len(nrow)), class="factor"
    )
    rowdf <- data.frame(
        ridx = seq_len(nrow),
        n = tabulate(x$ridx, nrow),
        sum = vapply(split(x$value, ridx), sum, numeric(1), USE.NAMES=FALSE)
    )

    ## columns: summarized cells (complete) in current sample
    ucidx <- unique(x$cidx)
    x$cidx <- match(x$cidx, ucidx)
    coldf <- data.frame(
        cidx = ucidx,
        n = tabulate(x$cidx, length(ucidx)),
        sum = vapply(split(x$value, x$cidx), sum, numeric(1), USE.NAMES=FALSE)
    )

    list(rowdf = rowdf, coldf = coldf)
}

The margin summary can be calculated as

register(MulticoreParam(progressbar=TRUE))
result <- tenxiterate(tenx, margin.summary, nrow = nrow(tenx))
rows <- Reduce(function(x, y) {
    idx <- c("n", "sum")
    x[, idx] <- x[, idx] + y[, idx]
    x
}, lapply(result, `[[`, 1))
cols <- do.call("rbind", lapply(result, `[[`, 2))

The summary takes about 8 minutes to read and process the entire million-cell data set using 6 cores and yieldSize = 10000.

Exploratory analysis

We return to our sampled SummarizedExperiment

se
table(se$Mouse)

With a reasonable subset of data in memory, it is possible to explore basic properties of the data.

The data is very sparse

sum(assay(se) == 0) / prod(dim(se))

Here are histograms of library size and reads per gene

hist(log10(1 + colSums(assay(se))))
hist(log(1 + rowSums(assay(se))))

Pooling across cells, the 'MA' plot is reassuringly familiar and approximately symmetric about Y = 0.

ma <- log(1 + rowsum(t(assay(se)), se$Mouse))
M <- ma[1,] - ma[2,]
A <- (ma[1,] + ma[2,]) / 2
plot(M ~ A)
abline(0, 0, lwd=2, col="blue")

Samples do not show obvious patterns with respect to mouse-of-origin.

library(Rtsne)
d <- dist(t(log(1 + assay(se))), method="manhattan")
tsne <- Rtsne(d)
plot(tsne$Y, pch=20, col = se$Mouse, cex=2, asp=1)

Session info

sessionInfo()