Running Banksy on Harmonized Seurat Object

[rstagnit]

Hi,

Banksy has worked well for me on individual samples and returned interesting results compared to traditional clustering. I am now interested in comparing the clustering results using Banksy on my harmonized data of 10 different samples to that I achieved with traditional seurat clustering. Following the tutorial (https://github.com/satijalab/seurat-wrappers/blob/master/docs/banksy.md#spatial-data-integration-with-harmony), I see that you start with a SpatialExperiment. I am wondering if it is possible to use a merged seurat object, similar to how I would with SCTransform. Using SCTransform, I would perform SCT on each sample individually, then merge the objects, then perform SCT again before running PCA and Harmony. However, I am having trouble accomplishing the same strategy with Banksy. When attempting to run Banksy on the merged object

DefaultAssay(Merged_postBanksy) <- "Spatial"
Merged_postBanksy <- RunBanksy(Merged_postBanksy, lambda = 0.2, verbose=TRUE, 
                               assay = 'Spatial', group = 'Sample', slot = 'data', features = 'all',
                               k_geom = 15)

I am met with the following error: Error in GetAssayData(): ! GetAssayData doesn't work for multiple layers in v5 assay.

This error is suggesting to use JoinLayers, but I did that to no success. I thought perhaps the issue is the lack of specified "dimx" or "dimy", but I wasn't sure if you had a suggestion.

Worst case, I can try doing this from a SpatialExperiment, but I am not familiar with using those objects and was similarly having difficulty.

Any assistance would be greatly appreciated.

Rob

[jleechung]

Hi @rstagnit,

Does calling JoinLayers after merging and before RunBanksy work? After merging your Seurat v5 objects, do you see something like this, where the counts and data from the objects are stored as separate layers?

> seu # merged seurat
An object of class Seurat 
120 features across 1500 samples within 1 assay 
Active assay: RNA (120 features, 0 variable features)
 6 layers present: counts.1.1, counts.2, counts.2.1, data.1.1, data.2.1, data.2

If so, running JoinLayers should merge the counts and data layers:

> JoinLayers(seu)
An object of class Seurat 
120 features across 1500 samples within 1 assay 
Active assay: RNA (120 features, 0 variable features)
 2 layers present: data, counts

Also, as you've pointed out, you'd need to supply the metadata column names corresponding to the spatial coordinates. It's best to use the same names (e.g. dimx, dimy) across all your Seurat objects before merging, so that calling merge row binds these columns together.

I'm not too familiar with the recommended SCTransform workflow. But following the SCT-Harmony workflow here, plugging in BANKSY may look something like:

# assumes your Seurat objects are stored in a list `seu_list`

# run SCTransform on individual Seurat objects in seu_list ...
seu = Reduce(merge, seu_list)
seu = JoinLayers(seu) # run this if the counts and data are stored in separate layers

seu = RunBanksy(
    seu, lambda = 0.2, assay = 'SCT', slot = 'data', group = 'Sample', 
    dimx = 'dimx', dimy = 'dimy',  k_geom = 15, split.scale = FALSE)

seu = RunPCA(seu, assay = 'BANKSY', features = rownames(seu))
seu = RunHarmony(seu, group.by.vars = 'Sample')
seu = RunUMAP(seu, reduction = 'harmony')
seu = FindNeighbors(seu, reduction = 'harmony')
seu = FindClusters(seu)

[rstagnit]

Hi @jleechung ,

I did as you suggested, running SCTransform on each individual object in seu_list then merged. The merged object is below.

> seu
An object of class Seurat 
36154 features across 32517 samples within 2 assays 
Active assay: SCT (18069 features, 0 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: Spatial
 2 images present: slice1, slice1.2
> seu = JoinLayers(seu)
Error in UseMethod(generic = "JoinLayers", object = object) : 
  no applicable method for 'JoinLayers' applied to an object of class "c('SCTAssay', 'Assay', 'KeyMixin')"

Running JoinLayers returns that error. Attempting to RunBanksy without using JoinLayers returns the following:

> seu = RunBanksy(
+   seu, lambda = 0.2, assay = 'SCT', slot = 'data', group = 'Sample', 
+   dimx = 'dimx', dimy = 'dimy',  k_geom = 15, split.scale = FALSE)
Fetching data from slot data from assay SCT
Subsetting by features
Error: SCT assay is comprised of multiple SCT models. To change the variable features, please set manually with VariableFeatures<-
In addition: Warning message:
In get_data(object, assay, slot, features, verbose) :
  No variable features found. Running Seurat::FindVariableFeatures

I then tried running the following:

> VariableFeatures(seu) <- c(VariableFeatures(AM_1), VariableFeatures(AM_2), VariableFeatures(AM_3), VariableFeatures(AM_6), VariableFeatures(AM_16), VariableFeatures(CM_1), VariableFeatures(CM_C10), VariableFeatures(CM_C20), VariableFeatures(CM_A26))
> seu = RunBanksy(
+   seu, lambda = 0.2, assay = 'SCT', slot = 'data', group = 'Sample', 
+   dimx = 'dimx', dimy = 'dimy',  k_geom = 15, split.scale = FALSE)
Fetching data from slot data from assay SCT
Subsetting by features
Staggering locations by Sample
Computing neighbors...
Spatial mode is kNN_median
Parameters: k_geom=15
Done
Computing harmonic m = 0
Using 15 neighbors
Done
Creating Banksy matrix
Scaling BANKSY matrix. Do not call ScaleData on assay BANKSY

Which my session got stuck here and ultimately caused R to abort.

[jleechung]

Thanks for the console output.

• It seems like JoinLayers doesn't need to be called, since counts and data are merged in the SCT assay.

• As for the stalled R session - could you try using a smaller number of features just to see if it's a memory / runtime issue? The Seurat team seems to recommend selecting 1) the intersection of the variable features of each dataset, or 2) by ranking features by the number of datasets in which they are deemed highly variable, and selecting features based on the ranks.

seu = Reduce(merge, seu_list) # merge

# option 1: intersection across samples
feat = Reduce(intersect, lapply(seu_list, VariableFeatures))
# option 2: top variable features across samples
feat = SelectIntegrationFeatures(seu_list, nfeatures = 2000)

VariableFeatures(seu) = feat # add variable features to merged seurat object

seu = RunBanksy(
    seu, lambda = 0.2, assay = 'SCT', slot = 'data', group = 'Sample', 
    dimx = 'dimx', dimy = 'dimy',  k_geom = 15, split.scale = FALSE)

[rstagnit]

Hi @jleechung ,

Sorry for the delayed reply. This method that you have recommended has worked, although now I was left with lots of clusters that only belong to one specific sample, which doesn't happen when I use traditional clustering. I have realized that SCtransform in its attempts to capture rare cell types, has a tendency to do that for integrated samples. I switch back to doing the traditional log normalize, find variable features, scale data approach with good results.

I want to point out a potential issue people may run into if they try to follow this exactly. When using Reduce(merge, seu_list), R has issues naming the slices as seen below.

> seu = Reduce(merge, seu_list) Warning: Some cell names are duplicated across objects provided. Renaming to enforce unique cell names. Warning: Some cell names are duplicated across objects provided. Renaming to enforce unique cell names. Warning: Some cell names are duplicated across objects provided. Renaming to enforce unique cell names. Warning: Some cell names are duplicated across objects provided. Renaming to enforce unique cell names. Warning: Some cell names are duplicated across objects provided. Renaming to enforce unique cell names. Warning messages: 1: Key ‘slice1_’ taken, using ‘slice12_’ instead 2: Key ‘slice1_’ taken, using ‘slice12_’ instead 3: Key ‘slice1_’ taken, using ‘slice12_’ instead 4: Key ‘slice1_’ taken, using ‘slice12_’ instead 5: Key ‘slice1_’ taken, using ‘slice12_’ instead 6: Key ‘slice1_’ taken, using ‘slice12_’ instead 7: Key ‘slice1_’ taken, using ‘slice12_’ instead 8: Key ‘slice1_’ taken, using ‘slice12_’ instead 9: Key ‘slice1_’ taken, using ‘slice12_’ instead 10: Key ‘slice1_’ taken, using ‘slice12_’ instead

This makes it so that only the first and last sample in the list have unique slice names, and if you try to run SpatialDimPlot, those will be the only two shown. You can circumvent this issue by just merging via: > test <-merge(M_1, y = c(M_2, M_3, M_6, M_16, M_C1, M_C10, M_C20, M_A26, M_C12, M_C18)) Warning: Some cell names are duplicated across objects provided. Renaming to enforce unique cell names. Warning messages: 1: Key ‘slice1_’ taken, using ‘slice12_’ instead 2: Key ‘slice1_’ taken, using ‘slice13_’ instead 3: Key ‘slice1_’ taken, using ‘slice14_’ instead 4: Key ‘slice1_’ taken, using ‘slice15_’ instead 5: Key ‘slice1_’ taken, using ‘slice16_’ instead 6: Key ‘slice1_’ taken, using ‘slice17_’ instead 7: Key ‘slice1_’ taken, using ‘slice18_’ instead 8: Key ‘slice1_’ taken, using ‘slice19_’ instead 9: Key ‘slice1_’ taken, using ‘slice110_’ instead 10: Key ‘slice1_’ taken, using ‘slice111_’ instead

As you can see, all samples now have unique names, and all will appear in the SpatialDimPlot.

Best, R

[jleechung]

Thanks @rstagnit, this is helpful and good to know! Closing this now, feel free to reopen or submit another issue if you run into any other problems.