modelGeneVar() extremely slow with HDF5-backed SingleCellExperiment

[keyingkuang]

Hi scran team,

Thanks for the great package! I’m encountering a performance issue when using modelGeneVar() on a large HDF5-backed SingleCellExperiment.

Previously, I used:
variable_genes <- scran::modelGeneVar(sce_object)
This completed in around 10 hours on the same large, HDF5-backed dataset. However, running the same call now takes several days without completing.

Here’s how we construct the object:
log_object <- HDF5Array::H5ADMatrix(raw_file_path, "data") count_object <- HDF5Array::H5ADMatrix(raw_file_path, "counts") sce_object <- SingleCellExperiment::SingleCellExperiment( assays = list(counts = count_object, logcounts = log_object) )

It seems that modelGeneVar() might not be fully optimized for delayed operations or is trying to load the entire matrix into memory. However, I don’t see an explicit memory error.

Do you have any suggestions on how to speed this up or safely apply modelGeneVar() to HDF5-backed data?

Thanks,
Maggie

[LTLA]

Oof. 10 hours already sounds pretty painful to me, never mind several days.

Anyway, short answer: consider using scrapper::modelGeneVariances() with beachmat.hdf5 installed.

library(scRNAseq)
sce <- fetchDataset("zilionis-lung-2019", "2023-12-20", path="human")
mat <- assay(sce)

sf <- scrapper::centerSizeFactors(colSums(mat))
# For simplicity, to avoid zero SFs from our un-QC'd data.
lmat <- scrapper::normalizeCounts(mat, size.factors=pmax(0.001, sf)) 
## DelayedArray::showtree(lmat)
## 41861x173954 double, sparse: DelayedMatrix object
## └─ 41861x173954 double, sparse: Set dimnames
##    └─ 41861x173954 double, sparse: Stack of 2 unary iso op(s)
##       └─ 41861x173954 double, sparse: Aperm (perm=c(2,1))
##          └─ 173954x41861 double, sparse: Unary iso op with args
##             └─ 173954x41861 integer, sparse: Aperm (perm=c(2,1))
##                └─ 41861x173954 integer, sparse: ReloadedArraySeed object
##                   └─ 41861x173954 integer, sparse: [seed] CSC_H5SparseMatrixSeed object

library(beachmat.hdf5)
system.time(res <- scrapper::modelGeneVariances(lmat, block=sce$Patient))
##    user  system elapsed
##   0.918   0.031   0.954

As for the long answer: the poor performance is because scran::modelGeneVar extracts blocks of rows for calculation of the row-wise variances. However, I'm assuming your HDF5 sparse matrix is stored in CSC format. (H5AD would declare it as CSR but that's because they store their genes in the columns, so from our perspective, it's CSC.) A CSC matrix is optimized for column access but not great for row access. For HDF5-backed matrices, this is even worse because each row access would have to read the entire file, over and over again, just to get each block of rows.

Now, why did it increase from 10 hours to several days+? I don't know. We just call various HDF5Array functions to extract the data, so I would guess that it's a change somewhere in that part of the ecosystem that is causing increased run-times. There were some quite major changes recently to the underlying DelayedArray + SparseArray packages, which may be the cause.

scrapper is a bit more optimized as it tries to use a more CSC-friendly access pattern when computing the variances (though this too has its inherent limits, see libscran/scrapper#10). scrapper was intended to be scran's high-performance replacement but we never got any funding to actually replace scran in Bioconductor documentation/workflows (e.g., OSCA book). So the package itself is complete but no one knows that they should be using it. :(