10x Genomics datasets as DelayedMatrix objects
TENxMatrix-class.RdA 10x Genomics dataset like the "1.3 Million Brain Cell Dataset" is an HDF5 sparse matrix stored in CSR/CSC/Yale format ("Compressed Sparse Row").
The TENxMatrix class is a DelayedMatrix subclass for representing and operating on this kind of dataset.
All the operations available for DelayedMatrix objects work on TENxMatrix objects.
Details
In addition to all the methods defined for
DelayedMatrix objects, TENxMatrix objects
support the following specialized methods: nzcount() and
extractNonzeroDataByCol(). See ?H5SparseMatrixSeed
for more information about what these methods do.
Note
If your dataset uses the HDF5 sparse matrix representation from
10x Genomics, use the TENxMatrix() constructor documented here.
But if your dataset uses the conventional (a.k.a. dense) HDF5
representation, use the HDF5Array() constructor instead.
See also
HDF5Array objects for representing conventional (a.k.a. dense) HDF5 datasets as DelayedArray objects.
DelayedMatrix objects in the DelayedArray package.
writeTENxMatrixfor writing a matrix-like object as an HDF5-based sparse matrix.The
TENxBrainDatadataset (in the TENxBrainData package).detectCoresfrom the parallel package.setAutoBPPARAMandsetAutoBlockSizein the DelayedArray package.colAutoGridandblockApplyin the DelayedArray package.The TENxMatrixSeed helper class.
h5lsto list the content of an HDF5 file.NumericList and IntegerList objects in the IRanges package.
SparseArray objects in the SparseArray package.
Examples
## ---------------------------------------------------------------------
## SIMPLE TENxMatrix EXAMPLE
## ---------------------------------------------------------------------
sm <- Matrix::rsparsematrix(10, 7, density=0.3)
M <- writeTENxMatrix(sm)
M
#> <10 x 7> sparse TENxMatrix object of type "double":
#> [,1] [,2] [,3] ... [,6] [,7]
#> [1,] 0.000 0.000 0.240 . 0.00 0.00
#> [2,] -0.016 0.000 0.000 . 0.00 0.00
#> [3,] 0.000 -0.940 0.540 . 0.00 -0.91
#> [4,] -1.900 0.000 1.900 . 0.00 0.00
#> [5,] -1.300 -0.097 0.000 . -0.83 0.18
#> [6,] -0.130 0.000 0.000 . 0.00 0.00
#> [7,] 0.000 0.000 0.000 . 0.00 0.00
#> [8,] 0.000 0.000 0.000 . 0.94 0.74
#> [9,] 0.000 0.000 0.000 . 0.00 0.00
#> [10,] 0.000 0.000 0.000 . 0.00 0.00
class(M) # TENxMatrix
#> [1] "TENxMatrix"
#> attr(,"package")
#> [1] "HDF5Array"
is(M, "DelayedMatrix") # TRUE
#> [1] TRUE
seed(M)
#> <10 x 7> sparse TENxMatrixSeed object of type "double":
#> # dirname: /tmp/RtmpYQEKxJ/HDF5Array_dump
#> # basename: autof48f0a8541.h5
#> # group: /HDF5ArrayAUTO00001
class(seed(M)) # TENxMatrixSeed
#> [1] "TENxMatrixSeed"
#> attr(,"package")
#> [1] "HDF5Array"
rhdf5::h5ls(path(M))
#> group name otype dclass dim
#> 0 / HDF5ArrayAUTO00001 H5I_GROUP
#> 1 /HDF5ArrayAUTO00001 data H5I_DATASET FLOAT 21
#> 2 /HDF5ArrayAUTO00001 indices H5I_DATASET INTEGER 21
#> 3 /HDF5ArrayAUTO00001 indptr H5I_DATASET INTEGER 8
#> 4 /HDF5ArrayAUTO00001 shape H5I_DATASET INTEGER 2
is_sparse(M) # TRUE
#> [1] TRUE
## Use coercion to load the full dataset into memory:
as.matrix(M) # as ordinary array (usually not recommended)
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.000 0.000 0.24 0.36 0.000 0.00 0.00
#> [2,] -0.016 0.000 0.00 0.00 0.000 0.00 0.00
#> [3,] 0.000 -0.940 0.54 0.00 -1.500 0.00 -0.91
#> [4,] -1.900 0.000 1.90 0.11 1.600 0.00 0.00
#> [5,] -1.300 -0.097 0.00 0.00 -0.053 -0.83 0.18
#> [6,] -0.130 0.000 0.00 0.00 0.000 0.00 0.00
#> [7,] 0.000 0.000 0.00 0.00 0.470 0.00 0.00
#> [8,] 0.000 0.000 0.00 0.00 -0.280 0.94 0.74
#> [9,] 0.000 0.000 0.00 0.00 0.000 0.00 0.00
#> [10,] 0.000 0.000 0.00 0.00 0.000 0.00 0.00
as(M, "dgCMatrix") # as dgCMatrix (brings back 'sm')
#> 10 x 7 sparse Matrix of class "dgCMatrix"
#>
#> [1,] . . 0.24 0.36 . . .
#> [2,] -0.016 . . . . . .
#> [3,] . -0.940 0.54 . -1.500 . -0.91
#> [4,] -1.900 . 1.90 0.11 1.600 . .
#> [5,] -1.300 -0.097 . . -0.053 -0.83 0.18
#> [6,] -0.130 . . . . . .
#> [7,] . . . . 0.470 . .
#> [8,] . . . . -0.280 0.94 0.74
#> [9,] . . . . . . .
#> [10,] . . . . . . .
as(M, "SparseArray") # as SparseArray object (most efficient)
#> <10 x 7 SparseMatrix> of type "double" [nzcount=21 (30%)]:
#> [,1] [,2] [,3] ... [,6] [,7]
#> [1,] 0.000 0.000 0.240 . 0.00 0.00
#> [2,] -0.016 0.000 0.000 . 0.00 0.00
#> [3,] 0.000 -0.940 0.540 . 0.00 -0.91
#> [4,] -1.900 0.000 1.900 . 0.00 0.00
#> [5,] -1.300 -0.097 0.000 . -0.83 0.18
#> [6,] -0.130 0.000 0.000 . 0.00 0.00
#> [7,] 0.000 0.000 0.000 . 0.00 0.00
#> [8,] 0.000 0.000 0.000 . 0.94 0.74
#> [9,] 0.000 0.000 0.000 . 0.00 0.00
#> [10,] 0.000 0.000 0.000 . 0.00 0.00
SparseArray(M) # equivalent to 'as(M, "SparseArray")'
#> <10 x 7 SparseMatrix> of type "double" [nzcount=21 (30%)]:
#> [,1] [,2] [,3] ... [,6] [,7]
#> [1,] 0.000 0.000 0.240 . 0.00 0.00
#> [2,] -0.016 0.000 0.000 . 0.00 0.00
#> [3,] 0.000 -0.940 0.540 . 0.00 -0.91
#> [4,] -1.900 0.000 1.900 . 0.00 0.00
#> [5,] -1.300 -0.097 0.000 . -0.83 0.18
#> [6,] -0.130 0.000 0.000 . 0.00 0.00
#> [7,] 0.000 0.000 0.000 . 0.00 0.00
#> [8,] 0.000 0.000 0.000 . 0.94 0.74
#> [9,] 0.000 0.000 0.000 . 0.00 0.00
#> [10,] 0.000 0.000 0.000 . 0.00 0.00
## ---------------------------------------------------------------------
## THE "1.3 Million Brain Cell Dataset" AS A DelayedMatrix OBJECT
## ---------------------------------------------------------------------
## The 1.3 Million Brain Cell Dataset from 10x Genomics is available
## via ExperimentHub:
library(ExperimentHub)
#> Loading required package: AnnotationHub
#> Loading required package: BiocFileCache
#> Loading required package: dbplyr
#>
#> Attaching package: ‘AnnotationHub’
#> The following object is masked from ‘package:Biobase’:
#>
#> cache
hub <- ExperimentHub()
query(hub, "TENxBrainData")
#> ExperimentHub with 8 records
#> # snapshotDate(): 2025-04-12
#> # $dataprovider: 10X Genomics
#> # $species: Mus musculus
#> # $rdataclass: character
#> # additional mcols(): taxonomyid, genome, description,
#> # coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
#> # rdatapath, sourceurl, sourcetype
#> # retrieve records with, e.g., 'object[["EH1039"]]'
#>
#> title
#> EH1039 | Brain scRNA-seq data, 'HDF5-based 10X Genomics' format
#> EH1040 | Brain scRNA-seq data, 'dense matrix' format
#> EH1041 | Brain scRNA-seq data, sample (column) annotation
#> EH1042 | Brain scRNA-seq data, gene (row) annotation
#> EH1689 | Brain scRNA-seq data 20k subset, 'HDF5-based 10x Genomics' format
#> EH1690 | Brain scRNA-seq data 20k subset, 'dense matrix' format
#> EH1691 | Brain scRNA-seq data 20k subset, sample (column) annotation
#> EH1692 | Brain scRNA-seq data 20k subset, gene (row) annotation
fname <- hub[["EH1039"]]
#> see ?TENxBrainData and browseVignettes('TENxBrainData') for documentation
#> loading from cache
## 'fname' is an HDF5 file. Use h5ls() to list its content:
h5ls(fname)
#> group name otype dclass dim
#> 0 / mm10 H5I_GROUP
#> 1 /mm10 barcodes H5I_DATASET STRING 1306127
#> 2 /mm10 data H5I_DATASET INTEGER 2624828308
#> 3 /mm10 gene_names H5I_DATASET STRING 27998
#> 4 /mm10 genes H5I_DATASET STRING 27998
#> 5 /mm10 indices H5I_DATASET INTEGER 2624828308
#> 6 /mm10 indptr H5I_DATASET INTEGER 1306128
#> 7 /mm10 shape H5I_DATASET INTEGER 2
## The 1.3 Million Brain Cell Dataset is represented by the "mm10"
## group. We point the TENxMatrix() constructor to this group to
## create a TENxMatrix object representing the dataset:
oneM <- TENxMatrix(fname, group="mm10")
oneM
#> <27998 x 1306127> sparse TENxMatrix object of type "integer":
#> AAACCTGAGATAGGAG-1 ... TTTGTCATCTGAAAGA-133
#> ENSMUSG00000051951 0 . 0
#> ENSMUSG00000089699 0 . 0
#> ENSMUSG00000102343 0 . 0
#> ENSMUSG00000025900 0 . 0
#> ENSMUSG00000109048 0 . 0
#> ... . . .
#> ENSMUSG00000079808 0 . 0
#> ENSMUSG00000095041 1 . 0
#> ENSMUSG00000063897 0 . 0
#> ENSMUSG00000096730 0 . 0
#> ENSMUSG00000095742 0 . 0
is(oneM, "DelayedMatrix") # TRUE
#> [1] TRUE
seed(oneM)
#> <27998 x 1306127> sparse TENxMatrixSeed object of type "integer":
#> # dirname: /github/home/.cache/R/ExperimentHub
#> # basename: 7b946c52a4b_1039
#> # group: /mm10
path(oneM)
#> [1] "/github/home/.cache/R/ExperimentHub/7b946c52a4b_1039"
nzcount(oneM) # nb of nonzero values in the dataset
#> [1] 2624828308
## Some examples of delayed operations:
oneM != 0
#> <27998 x 1306127> sparse DelayedMatrix object of type "logical":
#> AAACCTGAGATAGGAG-1 ... TTTGTCATCTGAAAGA-133
#> ENSMUSG00000051951 FALSE . FALSE
#> ENSMUSG00000089699 FALSE . FALSE
#> ENSMUSG00000102343 FALSE . FALSE
#> ENSMUSG00000025900 FALSE . FALSE
#> ENSMUSG00000109048 FALSE . FALSE
#> ... . . .
#> ENSMUSG00000079808 FALSE . FALSE
#> ENSMUSG00000095041 TRUE . FALSE
#> ENSMUSG00000063897 FALSE . FALSE
#> ENSMUSG00000096730 FALSE . FALSE
#> ENSMUSG00000095742 FALSE . FALSE
oneM^2
#> <27998 x 1306127> sparse DelayedMatrix object of type "double":
#> AAACCTGAGATAGGAG-1 ... TTTGTCATCTGAAAGA-133
#> ENSMUSG00000051951 0 . 0
#> ENSMUSG00000089699 0 . 0
#> ENSMUSG00000102343 0 . 0
#> ENSMUSG00000025900 0 . 0
#> ENSMUSG00000109048 0 . 0
#> ... . . .
#> ENSMUSG00000079808 0 . 0
#> ENSMUSG00000095041 1 . 0
#> ENSMUSG00000063897 0 . 0
#> ENSMUSG00000096730 0 . 0
#> ENSMUSG00000095742 0 . 0
## ---------------------------------------------------------------------
## SOME EXAMPLES OF ROW/COL SUMMARIZATION
## ---------------------------------------------------------------------
## In order to reduce computation times, we'll use only the first
## 25000 columns of the 1.3 Million Brain Cell Dataset:
oneM25k <- oneM[ , 1:25000]
## Row/col summarization methods like rowSums() use a block-processing
## mechanism behind the scene that can be controlled via global
## settings. 2 important settings that can have a strong impact on
## performance are the automatic number of workers and automatic block
## size, controlled by setAutoBPPARAM() and setAutoBlockSize()
## respectively.
library(BiocParallel)
if (.Platform$OS.type != "windows") {
## On a modern Linux laptop with 8 cores (as reported by
## parallel::detectCores()) and 16 Gb of RAM, reasonably good
## performance is achieved by setting the automatic number of workers
## to 5 or 6 and the automatic block size between 300 Mb and 400 Mb:
workers <- 5
block_size <- 3e8 # 300 Mb
setAutoBPPARAM(MulticoreParam(workers))
} else {
## MulticoreParam() is not supported on Windows so we use SnowParam()
## on this platform. Also we reduce the block size to 200 Mb on
## 32-bit Windows to avoid memory allocation problems (they tend to
## be common there because a process cannot use more than 3 Gb of
## memory).
workers <- 4
setAutoBPPARAM(SnowParam(workers))
block_size <- if (.Platform$r_arch == "i386") 2e8 else 3e8
}
setAutoBlockSize(block_size)
#> automatic block size set to 3e+08 bytes (was 1e+08)
## We're ready to compute the library sizes, number of genes expressed
## per cell, and average expression across cells:
system.time(lib_sizes <- colSums(oneM25k))
#> user system elapsed
#> 3.932 1.032 1.689
system.time(n_exprs <- colSums(oneM25k != 0))
#> user system elapsed
#> 11.985 1.559 4.333
system.time(ave_exprs <- rowMeans(oneM25k))
#> user system elapsed
#> 4.130 0.141 4.271
## Note that the 3 computations above load the data in oneM25k 3 times
## in memory. This can be avoided by computing the 3 summarizations in
## a single pass with blockApply(). First we define the function that
## we're going to apply to each block of data:
FUN <- function(block)
list(colSums(block), colSums(block != 0), rowSums(block))
## Then we call blockApply() to apply FUN() to each block. The blocks
## are defined by the grid passed to the 'grid' argument. In this case
## we supply a grid made with colAutoGrid() to generate blocks of full
## columns (see ?colAutoGrid for more information):
system.time({
block_results <- blockApply(oneM25k, FUN, grid=colAutoGrid(oneM25k),
verbose=TRUE)
})
#> user system elapsed
#> 8.345 2.376 7.638
## 'block_results' is a list with 1 list element per block in
## colAutoGrid(oneM25k). Each list element is the result that was
## obtained by applying FUN() on the block so is itself a list of
## length 3.
## Let's combine the results:
lib_sizes2 <- unlist(lapply(block_results, `[[`, 1L))
n_exprs2 <- unlist(lapply(block_results, `[[`, 2L))
block_rowsums <- unlist(lapply(block_results, `[[`, 3L), use.names=FALSE)
tot_exprs <- rowSums(matrix(block_rowsums, nrow=nrow(oneM25k)))
ave_exprs2 <- setNames(tot_exprs / ncol(oneM25k), rownames(oneM25k))
## Sanity checks:
stopifnot(all.equal(lib_sizes, lib_sizes2))
stopifnot(all.equal(n_exprs, n_exprs2))
stopifnot(all.equal(ave_exprs, ave_exprs2))
## Turn off parallel evaluation and reset automatic block size to factory
## settings:
setAutoBPPARAM()
setAutoBlockSize()
#> automatic block size set to 1e+08 bytes (was 3e+08)
## ---------------------------------------------------------------------
## extractNonzeroDataByCol()
## ---------------------------------------------------------------------
## extractNonzeroDataByCol() provides a convenient and very efficient
## way to extract the nonzero data in a compact form:
nonzeros <- extractNonzeroDataByCol(oneM, 1:25000) # takes < 5 sec.
## The data is returned as an IntegerList object with one list element
## per column and no row indices associated to the values in the object.
## Furthermore, the values within a given list element can be returned
## in any order:
nonzeros
#> IntegerList of length 25000
#> [[1]] 1 19 14 40 29 1 17 26 5 7 2 2 1 1 2 2 ... 5 1 1 2 2 1 1 1 1 1 1 14 1 1 1
#> [[2]] 6 1 7 16 15 10 13 5 1 1 1 1 1 1 1 1 1 ... 2 2 2 1 1 3 1 1 1 1 1 2 1 2 1 1
#> [[3]] 8 1 13 25 7 15 1 12 4 1 1 1 2 1 1 2 1 ... 2 3 2 1 2 1 1 3 1 1 2 7 2 1 1 1
#> [[4]] 2 1 11 1 6 2 18 14 12 8 4 2 1 1 1 1 3 ... 1 2 5 1 1 1 1 3 2 3 8 1 1 2 1 1
#> [[5]] 2 1 20 1 5 1 44 46 14 20 4 16 1 4 2 1 ... 1 2 2 1 1 3 1 3 20 1 1 2 2 1 2
#> [[6]] 1 2 41 1 4 33 2 73 102 44 55 12 34 2 1 ... 1 2 3 1 1 2 22 2 2 3 1 1 3 1 2
#> [[7]] 1 6 3 12 10 4 9 6 1 1 1 1 1 1 1 1 1 ... 1 1 6 1 1 1 1 2 1 2 1 1 1 6 1 1 2
#> [[8]] 1 1 20 8 23 27 1 15 13 3 4 1 1 3 1 1 ... 1 1 3 1 1 2 1 2 2 1 11 1 1 1 2 1
#> [[9]] 1 2 10 1 9 17 4 2 1 3 2 1 1 1 2 1 1 ... 1 1 1 2 1 1 1 2 1 2 1 1 1 8 1 1 1
#> [[10]] 3 18 14 2 1 29 28 1 14 23 4 15 3 1 2 1 ... 1 2 1 2 1 1 2 3 7 1 1 1 1 1 1
#> ...
#> <24990 more elements>
names(nonzeros) <- colnames(oneM25k)
## This can be used to compute some simple summaries like the library
## sizes and the number of genes expressed per cell. For these use
## cases, it is a lot more efficient than using colSums(oneM25k) and
## colSums(oneM25k != 0):
lib_sizes3 <- sum(nonzeros)
n_exprs3 <- lengths(nonzeros)
## Sanity checks:
stopifnot(all.equal(lib_sizes, lib_sizes3))
stopifnot(all.equal(n_exprs, n_exprs3))