Recommendations for very large SCE objects

[Thapeachydude]

Hello,

I hoping for some recommendations on how to handle very large SCE objects. I've seen the Chapter 14 Dealing with big data in the OSCA book, but I'm encountering an issue not mentioned there.

Specifically, I'd like to combine several SCEs but the size exeed the limit of 2^31-1 of an R matrix.
Converting the counts to a DelayedMatrix seems to circumvent the problem, I was just wondering what the "best" format would be....
HDF5SummarizedExperiment, TENxMatrix or another format?

I've also seen the recommendation of SVT_SparseMatrix in a recent issue, but many SCE functions such as logNormCounts or modelGeneVarseem to throw an error when dealing with an SVT_SparseMatrix.

Any insights would be much appreciated, thanks!

[LTLA]

If you can keep it in memory, I would use SVT_SparseMatrix. Newer packages like scrapper should handle these without problems. I would have expected that scran and scater can still take SVT_SparseMatrix objects without erroring out, so file an issue if that's not the case. (That said, the older packages are probably not efficient at processing them, I haven't gotten around to optimizing them.)

If you can't afford to keep it in memory, I would recommend using TileDBArray, which provides a TileDB-backed store that directly handles sparse data. TileDB enables efficient row- and column-based access to the same matrix. This is more useful than the compressed sparse formats we use with HDF5 (e.g., TENxMatrix), where CSC and CSR matrices can't be efficiently accessed by row and column, respectively.

[Thapeachydude]

Hi, thanks for the quick reply!
I can keep them in memory (in fact this is preferred, since I have to move them around a bit between servers and HDF5 always cause issues). I'll try the SVT_SparseMatrix + scrapper then. Is there any chance that you can provide a Bioconductor 3.18 compatible version?

[LTLA]

I suspect that you would need to do a lot of work to get the latest versions of these packages to work on 3.18. You could try, but there's a lot of dependencies that would need manual installation (e.g., assorthead, beachmat, BiocNeighbors, and then scrapper). The installation of the latest BiocNeighbors would probably break a few other packages in 3.18, ditto for beachmat.

In fact, I'm not even sure the SVT_SparseMatrix in 3.18 has the same definition as the class of the same name in 3.20+. IIRC there was a breaking change around that time, though I don't quite remember - something to do with the indices/values being swapped at the leaf nodes. You might need to manually install the latest SparseArray, DelayedArray and their dependencies as well.

IMO probably easiest to just switch to 3.20+.

[Thapeachydude]

I've upgraded to Bioconductor 3.20 and was able to convert the counts to an SVT_SparseMatrix.

<5001 x 8852792 SparseMatrix> of type "double" [nzcount=2894467569 (6.5%)]:
                  aaaaadia-1   aaaaajnp-1   aaaabcjp-1 ... oilknmkl-1 oilkooch-1
ENSG00000166535                     0                     0                     0   .                     0                     0
ENSG00000127837                     0                     0                     0   .                     0                     0
ENSG00000131043                     0                     0                     0   .                     0                     0
ENSG00000266967                     0                     0                     0   .                     0                     0
ENSG00000183044                     0                     0                     0   .                     0                     0
            ...                     .                     .                     .   .                     .                     .
ENSG00000197062                     0                     0                     0   .                     0                     0
ENSG00000130449                     0                     0                     0   .                     0                     0
ENSG00000153975                     0                     0                     0   .                     0                     0
ENSG00000162378                     0                     0                     0   .                     0                     0
ENSG00000159840                     1                     0                     0   .                     0                     0

After getting the size factors I run the normalization.
res <- normalizeCounts(counts(sce), size.factors = sizeFactors(sce))

This returns a DelayedArray. Trying to conver it to an SVT_SparseMatrix throws an error.

logcounts(sce) <- as(res, "SVT_SparseMatrix")
Error in .Call2(.NAME, ..., PACKAGE = "SparseArray") : 
  SVT_SparseArray object contains too many nonzero values to be turned into a COO_SparseArray object

Is this conversion of the logcounts necessary or is the DelayedArray fine to stay as is? (It doesn't seem to be written to a file).

[LTLA]

You can leave it as it is, scrapper doesn't care.

Seems like a bug in the coercion method, though, as normalizeCounts() doesn't change the number of non-zero elements. I'd suggest writing up an MRE. For example, this works for me:

library(SparseArray)
y <- as(abs(Matrix::rsparsematrix(10000, 1000, 0.1) * 10), "SVT_SparseMatrix")
y2 <- scrapper::normalizeCounts(y, size.factors=centerSizeFactors(colSums(y)))
class(as(y2, "SVT_SparseMatrix"))
## [1] "SVT_SparseMatrix"
## attr(,"package")
## [1] "SparseArray"

[Thapeachydude]

Apologies for the late reply. The bug persists.
I've tried to make an example with the code snippet above, but haven't been able to reproduce it that way.
The fact that it also disappears on smaller subsets of the "real" dataset is all the more annoying.

[LTLA]

Yeah I guess the error message:

  SVT_SparseArray object contains too many nonzero values to be turned into a COO_SparseArray object

says it all, now that I look at it. There is a COO_SparseArray intermediate that can't support the size of the matrix, but is otherwise unnecessary for the conversion from a dgCMatrix into a SVT_SparseMatrix. I dunno whether this is still present in 3.21, pinging @hpages for his thoughts.

[hpages]

Many coercion methods to SVT_SparseMatrix have been refactored in BioC 3.21, including coercion from DelayedMatrix to SVT_SparseMatrix. Most of them now avoid the intermediate COO_SparseArray representation so I would expect that upgrading to BioC 3.21 would solve the issue.

[hpages]

But if we are talking about in-memory representation, using a DelayedArray object doesn't really help here. It's cool that scrapper::normalizeCounts() returns instantly by delaying everything, but sooner or later you have to pay the bill. In this case, you pay it when you coerce back to SVT_SparseMatrix, and the price is much higher than with a normalizeCounts() that operates natively on the input matrix:

library(SparseArray)
library(SingleCellExperiment)
library(scrapper)

x <- poissonSparseMatrix(155000, 22000, density=0.1)

class(x)
# [1] "SVT_SparseMatrix"
# attr(,"package")
# [1] "SparseArray"

sf <- centerSizeFactors(colSums(x))

Going thru the DelayedArray representation:

system.time(normed0 <- scrapper::normalizeCounts(x, sf))  # cool, everything is delayed
#    user  system elapsed 
#   0.019   0.002   0.021 

class(normed0)
# [1] "DelayedMatrix"
# attr(,"package")
# [1] "DelayedArray"

system.time(normed1 <- as(normed0, "SVT_SparseMatrix"))  # time to pay the bill!
#    user  system elapsed 
#  50.671   7.331  58.042 

Operating natively on the SVT_SparseMatrix object:

native_normalizeCounts <- function(x, size.factors, log.base=2)
{   
    normalized <- t(t(x) / size.factors)
    log1p(normalized) / log(log.base)
}

system.time(normed2 <- native_normalizeCounts(x, sf))
#    user  system elapsed 
#  20.626   3.434  29.203 

identical(normed1, normed2)
# [1] TRUE

Note that native_normalizeCounts() will become even more efficient once the double transposition trick can be replaced with x / as_tile(size.factors, along=2). This is not ready yet but can already be emulated with SparseArray:::.Arith_SVT1_v2("/", x, size.factors, recycle.along=2L):

fast_native_normalizeCounts <- function(x, size.factors, log.base=2)
{   
    normalized <- SparseArray:::.Arith_SVT1_v2("/", x, size.factors, recycle.along=2L)
    log1p(normalized) / log(log.base)
}

system.time(normed3 <- fast_native_normalizeCounts(x, sf))
#    user  system elapsed 
#   4.725   0.173   5.698 

identical(normed1, normed3)
# [1] TRUE

fast_native_normalizeCounts() also significantly reduces memory usage.

x / as_tile(size.factors, along=2) will use SparseArray:::.Arith_SVT1_v2("/", x, size.factors, recycle.along=2L) behind the scene and should be ready soon in BioC 3.21, hopefully in the next couple of weeks.

H.

[hpages]

@vjcitn

[Thapeachydude]

As a quick follow-up: I finally got around to building a singularity image with R4.5.1 and BioC 3.21.
Unfortunately, converting the logcounts into an SVT_SparseMatrix still throws an error:

Error in .Call2(.NAME, ..., PACKAGE = "SparseArray") :
  build_OPBufTree_from_Lindex2() not ready yet
Calls: lapply ... .subassign_SVT_by_Lindex -> SparseArray.Call -> .Call2

[hpages]

Should be addressed in SparseArray 1.9.1 (BioC 3.22). See Bioconductor/SparseArray@0a2c62a

[hpages]

Also pushed to SparseArray 1.8.1 (BioC 3.21). See Bioconductor/SparseArray@ca78384

Both SparseArray 1.8.1 and 1.9.1 should become available via BiocManager::install() in the next couple of days.

[Thapeachydude]

Hi, I've updated to SparseArray 1.8.1 and now get:

 *** caught segfault ***
address 0x14eeec7cf000, cause 'memory not mapped'
Error in .Call2(.NAME, ..., PACKAGE = "SparseArray") :
  long vectors not supported yet: ../../src/include/Rinlinedfuns.h:551
Calls: lapply ... .subassign_SVT_by_Lindex -> SparseArray.Call -> .Call2
Execution halted