Memory and Time Issues with Xenium 5K samples

[atrull314]

Hi all,

We’ve been using BANKSY to perform analysis on Xenium samples. We’ve so far used it on two Xenium projects both of which used the standard ~300-400 gene panel, however, we’ve recently begun working with a project that contains 12 samples using the new 5K panel, and we’ve noticed a major increase in time and memory usage.

As background, we’ve normalized and integrated these samples and are running these through BANKSY with a number of combinations of lambdas, k values, pc values, and resolutions so we can determine the best values for this dataset. We’ve separated out the various BANKSY function calls into separate scripts to allow us to parallelize on the HPC we execute analyses on and to help determine where any memory/time bottlenecks are. We have found that executing clusterBanksy on this dataset for a single parameter combination has a peak memory at 512GB and the time to run this function can peak at 53 hours. This is causing a noticeable bottleneck in our workflow as a result, so we wanted to ask if there any plans on further optimization or if there were any additional suggestions on what we could do?

As some additional information, we are currently using the following 'flow' in order to execute Banksy:

computeBanksy() > runBanksyPCA() > RunHarmony() > runBanksyUMAP() > clusterBanksy()

And here's the code we are using to execute clusterBanksy (I'm happy to provide additional code snippets and any logging that may be informative). I'll also add that we are using SpatialExperiment objects to execute Banksy, and these have been converted from Seurat objects prior to starting the Banksy portion of the analysis.

spe_xenium_obj <- clusterBanksy(
    spe_xenium_obj,
    dimred = "BANKSY_harmony",
    resolution = opt$res,
    lambda = opt$lambda,
    ndims = opt$nPCs,
    seed = 1234,
    verbose = TRUE
)

[jleechung]

Hi @atrull314,

The easiest thing worth trying is using the SeuratWrappers interface for Banksy. This leverages all of Seurat's optimizations which should speed things up. You wouldn't need to convert your Seurat objects to SpatialExperiment objects as well.

You can do a parameter sweep that parallelizes over lambda and k_geom, and within each job, sweep over number of PCs and cluster resolution. Here is a sketch:

# Here seu is a combined Seurat objects comprising all samples, 
# with some metadata column `sample_id` denoting the sample id
# Parallelize across lambda, k_geom combinations
seu <- RunBanksy(seu, 
                 lambda = opt$lambda,
                 k_geom = opt$k_geom,
                 assay = "Spatial",
                 slot = "data",
                 features = "variable",
                 group = "sample_id",
                 split.scale = TRUE)
seu <- RunPCA(seu, assay = "BANKSY", npcs = 50)
# Run Harmony on all PCs, then sweep over nPCs
# One could also run Harmony for different nPCs
seu <- RunHarmony(seu, group.by.vars = "sample_id", reduction = "pca")
for (npc in c(20, 30, 40)) {
  seu <- FindNeighbors(seu, reduction = "harmony", dims = 1:npc)
  seu <- FindClusters(seu, resolution = c(0.4, 0.8, 1.2, 1.6),
                      cluster.name = paste0("clusters_pc", npc))
}

Check out the link above, happy to answer any questions you might have.

[atrull314]

Awesome, thanks for the help! We're currently testing this and will let you know if anything comes up!

[atrull314]

Hi all,

Just wanted to follow up. Thank you for the suggestion -- it does look like that the SeuratWrapper has been able to significantly decrease runtime with our initial testing.

However, we have noted that the RunBanksy() function is still peaking at around 500GB of memory with the Xenium 5K data. I'll note that prior to running RunBanksy(), the Seurat object is filtered to the top 2000 variable features, so the 500GB peak occurs even after the filtering, which is also true of the workflow I mentioned in the first reply.

And with regards to the initial post, I did also go through the clusterBanksy() code and saw the use of the leidenAlg package for the clustering algorithm. I'm not sure if it might be worth trying to swap that for a package that supported threading to help with time?

[jleechung]

How many cells do you have across the 12 samples? The peak would most certainly correspond to the memory required to store the scaled BANKSY matrix in its full size. We have had several discussions on memory handling for large datasets (see #39 for instance), but none of this has been implemented within the official Banksy or SeuratWrappers codebases yet.

As for the use of leidenAlg -- it's certainly worth trying to swap that, though we do not have bandwidth to implement that at the moment. You're welcome to raise a PR!