DelayedMatrixStats ports the matrixStats API to work with DelayedMatrix objects from the DelayedArray package. It provides high-performing functions operating on rows and columns of DelayedMatrix objects, including all subclasses such as RleArray (from the DelayedArray package) and HDF5Array (from the HDF5Array) as well as supporting all types of seeds, such as matrix (from the base package) and Matrix (from the Matrix package).
The DelayedArray package allows developers to store array-like data using in-memory or on-disk representations (e.g., in HDF5 files) and provides a common and familiar array-like interface for interacting with these data.
The DelayedMatrixStats package is designed to make life easier for Bioconductor developers wanting to use DelayedArray by providing a rich set of column-wise and row-wise summary functions.
We briefly demonstrate and explain these two features using a simple example. We’ll simulate some (unrealistic) RNA-seq read counts data from 10,000 genes and 20 samples and store it on disk as a HDF5Array:
library(DelayedArray)
x <- do.call(cbind, lapply(1:20, function(j) {
rpois(n = 10000, lambda = sample(20:40, 10000, replace = TRUE))
}))
colnames(x) <- paste0("S", 1:20)
x <- realize(x, "HDF5Array")
x
#> <10000 x 20> matrix of class HDF5Matrix and type "integer":
#> S1 S2 S3 S4 ... S17 S18 S19 S20
#> [1,] 24 37 48 31 . 26 33 17 33
#> [2,] 24 39 33 44 . 28 27 32 22
#> [3,] 21 21 30 15 . 27 26 18 23
#> [4,] 41 33 17 38 . 36 26 33 29
#> [5,] 24 27 34 20 . 25 35 26 37
#> ... . . . . . . . . .
#> [9996,] 27 37 21 37 . 19 18 22 38
#> [9997,] 40 44 38 21 . 20 21 25 36
#> [9998,] 40 27 37 48 . 36 24 23 20
#> [9999,] 29 36 26 25 . 37 30 47 36
#> [10000,] 23 21 19 33 . 31 21 40 23Suppose you wish to compute the standard deviation of the read counts for each gene.
You might think to use apply() like in the following:
system.time(row_sds <- apply(x, 1, sd))
#> user system elapsed
#> 128.88 15.86 144.74
head(row_sds)
#> [1] 8.125820 7.720615 9.414994 7.265455 8.443465 8.642003This works, but takes quite a while.
Or perhaps you already know that the matrixStats package
provides a rowSds() function:
matrixStats::rowSds(x)
#> Error in rowVars(x, rows = rows, cols = cols, na.rm = na.rm, center = center, : Argument 'x' must be a matrix or a vector.Unfortunately (and perhaps unsurprisingly) this doesn’t work. matrixStats is designed for use on in-memory matrix objects. Well, why don’t we just first realize our data in-memory and then use matrixStats
system.time(row_sds <- matrixStats::rowSds(as.matrix(x)))
#> user system elapsed
#> 0.01 0.00 0.01
head(row_sds)
#> [1] 8.125820 7.720615 9.414994 7.265455 8.443465 8.642003This works and is many times faster than the apply()-based approach! However,
it rather defeats the purpose of using a HDF5Array for storing the
data since we have to bring all the data into memory at once to compute the
result.
Instead, we can use DelayedMatrixStats::rowSds(), which has the speed
benefits of matrixStats::rowSds()1 In fact, it currently uses matrixStats::rowSds() under the hood. but without having to load the
entire data into memory at once2 In this case, it loads blocks of data row-by-row. The amount of
data loaded into memory at any one time is controlled by the
default block size global setting; see ?DelayedArray::getAutoBlockSize
for details. Notably, if the data are small enough (and the default block size
is large enough) then all the data is loaded as a single block, but this
approach generalizes and still works when the data are too large to be
loaded into memory in one block.:
library(DelayedMatrixStats)
system.time(row_sds <- rowSds(x))
#> user system elapsed
#> 0.01 0.00 0.02
head(row_sds)
#> [1] 8.125820 7.720615 9.414994 7.265455 8.443465 8.642003Finally, by using DelayedMatrixStats we can use the same code,
(colMedians(x)) regardless of whether the input is an ordinary matrix or a
DelayedMatrix. This is useful for packages wishing to support both types of
objects, e.g., packages wanting to retain backward compatibility or during a
transition period from matrix-based to DelayeMatrix-based objects.
The initial release of DelayedMatrixStats supports the complete column-wise and row-wise API matrixStats API3 Some of the API is covered via inheritance to functionality in DelayedArray. Please see the matrixStats vignette (available online) for a summary these methods. The following table documents the API coverage and availability of ‘seed-aware’ methods in the current version of DelayedMatrixStats, where:
| Method | Block processing | base::matrix optimized | Matrix::dgCMatrix optimized | Matrix::lgCMatrix optimized | DelayedArray::RleArray (SolidRleArraySeed) optimized | DelayedArray::RleArray (ChunkedRleArraySeed) optimized | HDF5Array::HDF5Matrix optimized | base::data.frame optimized | S4Vectors::DataFrame optimized |
|---|---|---|---|---|---|---|---|---|---|
colAlls() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colAnyMissings() |
<U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colAnyNAs() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colAnys() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colAvgsPerRowSet() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colCollapse() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colCounts() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colCummaxs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colCummins() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colCumprods() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colCumsums() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colIQRDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colIQRs() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colLogSumExps() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colMadDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colMads() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colMaxs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colMeans2() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colMedians() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colMins() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colOrderStats() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colProds() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colQuantiles() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colRanges() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colRanks() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colSdDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colSds() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colSums2() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colTabulates() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colVarDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colVars() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colWeightedMads() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colWeightedMeans() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colWeightedMedians() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colWeightedSds() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colWeightedVars() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
colsum() |
<U+2611><U+FE0F> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowAlls() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowAnyMissings() |
<U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowAnyNAs() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowAnys() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowAvgsPerColSet() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowCollapse() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowCounts() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowCummaxs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowCummins() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowCumprods() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowCumsums() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowIQRDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowIQRs() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowLogSumExps() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowMadDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowMads() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowMaxs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowMeans2() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowMedians() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowMins() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowOrderStats() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowProds() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowQuantiles() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowRanges() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowRanks() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowSdDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowSds() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowSums2() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowTabulates() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowVarDiffs() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowVars() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowWeightedMads() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowWeightedMeans() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowWeightedMedians() |
<U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowWeightedSds() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowWeightedVars() |
<U+2714> | <U+2714> | <U+2714> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
rowsum() |
<U+2611><U+FE0F> | <U+274C> | <U+2714> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> | <U+274C> |
As well as offering a familiar API, DelayedMatrixStats provides ‘seed-aware’ methods that are optimized for specific types of DelayedMatrix objects.
To illustrate this idea, we will compare two ways of computing the column sums of a DelayedMatrix object:
force_block_processing argumentWe will demonstrate this by computing the column sums matrices with 20,000 rows and 600 columns where the data have different structure and are stored in DelayedMatrix objects with different types of seed:
We use the microbenchmark package to measure running time and the profmem package to measure the total memory allocations of each method.
In each case, the ‘seed-aware’ method is many times faster and allocates substantially lower total memory.
library(DelayedMatrixStats)
library(sparseMatrixStats)
library(microbenchmark)
library(profmem)
set.seed(666)
# -----------------------------------------------------------------------------
# Dense with values in (0, 1)
# Fast, memory-efficient column sums of DelayedMatrix with ordinary matrix seed
#
# Generate some data
dense_matrix <- matrix(runif(20000 * 600),
nrow = 20000,
ncol = 600)
# Benchmark
dm_matrix <- DelayedArray(dense_matrix)
class(seed(dm_matrix))
#> [1] "matrix" "array"
dm_matrix
#> <20000 x 600> matrix of class DelayedMatrix and type "double":
#> [,1] [,2] [,3] ... [,599] [,600]
#> [1,] 0.7743685 0.6601787 0.4098798 . 0.89118118 0.05776471
#> [2,] 0.1972242 0.8436035 0.9198450 . 0.31799523 0.63099417
#> [3,] 0.9780138 0.2017589 0.4696158 . 0.31783791 0.02830454
#> [4,] 0.2013274 0.8797239 0.6474768 . 0.55217184 0.09678816
#> [5,] 0.3612444 0.8158778 0.5928599 . 0.08530977 0.39224147
#> ... . . . . . .
#> [19996,] 0.19490291 0.07763570 0.56391725 . 0.09703424 0.62659353
#> [19997,] 0.61182993 0.01910121 0.04046034 . 0.59708388 0.88389731
#> [19998,] 0.12932744 0.21155070 0.19344085 . 0.51682032 0.13378223
#> [19999,] 0.18985573 0.41716539 0.35110782 . 0.62939661 0.94601427
#> [20000,] 0.87889047 0.25308041 0.54666920 . 0.81630322 0.73272217
microbenchmark(
block_processing = colSums2(dm_matrix, force_block_processing = TRUE),
seed_aware = colSums2(dm_matrix),
times = 10)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> block_processing 71.23879 76.60649 155.66325 87.16033 95.79196 451.33492 10
#> seed_aware 16.96937 17.14279 17.31618 17.33544 17.49904 17.66521 10
#> cld
#> b
#> a
total(profmem(colSums2(dm_matrix, force_block_processing = TRUE)))
#> [1] 96260576
total(profmem(colSums2(dm_matrix)))
#> [1] 166120
# -----------------------------------------------------------------------------
# Sparse (60% zero) with values in (0, 1)
# Fast, memory-efficient column sums of DelayedMatrix with ordinary matrix seed
#
# Generate some data
sparse_matrix <- dense_matrix
zero_idx <- sample(length(sparse_matrix), 0.6 * length(sparse_matrix))
sparse_matrix[zero_idx] <- 0
# Benchmark
dm_dgCMatrix <- DelayedArray(Matrix(sparse_matrix, sparse = TRUE))
class(seed(dm_dgCMatrix))
#> [1] "dgCMatrix"
#> attr(,"package")
#> [1] "Matrix"
dm_dgCMatrix
#> <20000 x 600> sparse matrix of class DelayedMatrix and type "double":
#> [,1] [,2] [,3] ... [,599] [,600]
#> [1,] 0.7743685 0.0000000 0.0000000 . 0.89118118 0.00000000
#> [2,] 0.1972242 0.0000000 0.9198450 . 0.00000000 0.00000000
#> [3,] 0.9780138 0.0000000 0.4696158 . 0.31783791 0.00000000
#> [4,] 0.0000000 0.8797239 0.6474768 . 0.55217184 0.00000000
#> [5,] 0.3612444 0.0000000 0.0000000 . 0.08530977 0.39224147
#> ... . . . . . .
#> [19996,] 0.1949029 0.0776357 0.0000000 . 0.09703424 0.00000000
#> [19997,] 0.0000000 0.0000000 0.0000000 . 0.00000000 0.88389731
#> [19998,] 0.0000000 0.2115507 0.1934408 . 0.00000000 0.00000000
#> [19999,] 0.1898557 0.0000000 0.3511078 . 0.62939661 0.94601427
#> [20000,] 0.8788905 0.2530804 0.0000000 . 0.00000000 0.73272217
microbenchmark(
block_processing = colSums2(dm_dgCMatrix, force_block_processing = TRUE),
seed_aware = colSums2(dm_dgCMatrix),
times = 10)
#> Unit: milliseconds
#> expr min lq mean median uq
#> block_processing 297.253958 309.562153 428.716401 315.275634 665.832993
#> seed_aware 7.631218 7.718786 7.871794 7.871922 7.986615
#> max neval cld
#> 749.086945 10 b
#> 8.249745 10 a
total(profmem(colSums2(dm_dgCMatrix, force_block_processing = TRUE)))
#> [1] 249644616
total(profmem(colSums2(dm_dgCMatrix)))
#> [1] 7400
# -----------------------------------------------------------------------------
# Dense with values in {0, 100} featuring runs of identical values
# Fast, memory-efficient column sums of DelayedMatrix with Rle-based seed
#
# Generate some data
runs <- rep(sample(100, 500000, replace = TRUE), rpois(500000, 100))
runs <- runs[seq_len(20000 * 600)]
runs_matrix <- matrix(runs,
nrow = 20000,
ncol = 600)
# Benchmark
dm_rle <- RleArray(Rle(runs),
dim = c(20000, 600))
class(seed(dm_rle))
#> [1] "SolidRleArraySeed"
#> attr(,"package")
#> [1] "DelayedArray"
dm_rle
#> <20000 x 600> matrix of class RleMatrix and type "integer":
#> [,1] [,2] [,3] [,4] ... [,597] [,598] [,599] [,600]
#> [1,] 20 65 9 50 . 43 74 85 14
#> [2,] 20 65 9 50 . 43 74 85 14
#> [3,] 20 65 9 50 . 43 74 85 14
#> [4,] 20 65 9 50 . 43 74 85 14
#> [5,] 20 65 9 50 . 43 74 85 14
#> ... . . . . . . . . .
#> [19996,] 65 9 50 16 . 74 85 14 55
#> [19997,] 65 9 50 16 . 74 85 14 55
#> [19998,] 65 9 50 16 . 74 85 14 55
#> [19999,] 65 9 50 16 . 74 85 14 55
#> [20000,] 65 9 50 16 . 74 85 14 55
microbenchmark(
block_processing = colSums2(dm_rle, force_block_processing = TRUE),
seed_aware = colSums2(dm_rle),
times = 10)
#> Unit: milliseconds
#> expr min lq mean median uq
#> block_processing 1068.44901 1071.258866 1091.648249 1076.917459 1116.377664
#> seed_aware 6.25448 6.654302 9.134351 7.075909 7.536388
#> max neval cld
#> 1142.76207 10 b
#> 28.56634 10 a
total(profmem(colSums2(dm_rle, force_block_processing = TRUE)))
#> [1] 144416528
total(profmem(colSums2(dm_rle)))
#> [1] 35872The development of ‘seed-aware’ methods is ongoing work (see the Roadmap), and for now only a few methods and seed-types have a ‘seed-aware’ method.
An extensive set of benchmarks is under development at http://peterhickey.org/BenchmarkingDelayedMatrixStats/.
A key feature of a DelayedArray is the ability to register ‘delayed
operations’. For example, let’s compute sin(dm_matrix):
system.time(sin_dm_matrix <- sin(dm_matrix))
#> user system elapsed
#> 0 0 0This instantaneous because the operation is not actually performed, rather
it is registered and only performed when the object is realized. All methods
in DelayedMatrixStats will correctly realise these delayed
operations before computing the final result. For example, let’s compute
colSums2(sin_dm_matrix) and compare check we get the correct answer:
all.equal(colSums2(sin_dm_matrix), colSums(sin(as.matrix(dm_matrix))))
#> [1] TRUEThe initial version of DelayedMatrixStats provides complete coverage of the matrixStats column-wise and row-wise API4 Some of the API is covered via inheritance to functionality in DelayedArray, allowing package developers to use these functions with DelayedMatrix objects as well as with ordinary matrix objects. This should simplify package development and assist authors to support to their software for large datasets stored in disk-backed data structures such as HDF5Array. Such large datasets are increasingly common with the rise of single-cell genomics.
Future releases of DelayedMatrixStats will improve the
performance of these methods, specifically by developing additional ‘seed-aware’
methods. The plan is to prioritise commonly used methods (e.g.,
colMeans2()/rowMeans2(), colSums2()/rowSums2(), etc.) and the
development of ‘seed-aware’ methods for the HDF5Matrix class. To do so, we
will leverage the beachmat package. Proof-of-concept code
has shown that this can greatly increase the performance when analysing such
disk-backed data.
Importantly, all package developers using methods from DelayedMatrixStats will immediately gain from performance improvements to these low-level routines. By using DelayedMatrixStats, package developers will be able to focus on higher level programming tasks and address important scientific questions and technological challenges in high-throughput biology.