TENxGenomics.Rmd

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

```{r vignette_setup, echo=FALSE}
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.

```{r setup1, eval=FALSE}
biocLite("mtmorgan/TENxGenomics")
library(TENxGenomics)
```

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

```{r bfc}
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.

```{r}
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.

```{r}
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_.

```{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.

```{r}
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.

```{r}
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".

```{r}
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.

```{r}
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

```{r bfc-2}
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.

```{r}
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.

```{r}
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

```{r, eval=FALSE}
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

```{r}
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

```{r}
sum(assay(se) == 0) / prod(dim(se))
```

Here are histograms of library size and reads per gene

```{r}
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.

```{r}
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.

```{r}
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

```{r}
sessionInfo()
```

[TENxGenomics]: https:/github.com/mtmorgan/TENxGenomics
[10xGenomics]: https://support.10xgenomics.com/single-cell/datasets
[BiocFileCache]: https://bioconductor.org/packages/BiocFileCache
[SummarizedExperiment]: https://bioconductor.org/packages/SummarizedExperiment
[DelayedArray]: https://bioconductor.org/packages/DelayedArray