README.md
In mtmorgan/TENxGenomics: Interface to 10x Genomics data
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.
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
```{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.
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
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()
mtmorgan/TENxGenomics documentation built on May 23, 2019, 8:19 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.
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
```{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.
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
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()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.