Question after running FindAllMarkers on banksy clusters in CosMx data

[Chiranjit1504]

Hi,

Firstly, thank you for developing the BANKSY package — it is extremely user-friendly and integrates very well with Seurat.

I am currently working with a spatial dataset comprising multiple samples across several cohorts of liver diseases, and I am using BANKSY to identify spatial neighbourhoods. I have followed the suggested workflow from the vignette (BANKSY → PCA → Harmony; example attached below) and, to emphasise neighbourhood-level structure, I have set the hyperparameter λ = 0.8.

As a next step, I am interested in identifying genes that are differentially expressed across these neighbourhoods. For this, I ran FindAllMarkers using the “BANKSY” assay. While I obtain a list of marker genes, I notice that each gene is annotated with a suffix such as “m0” or “m1”. In some clusters, both suffixes appear, whereas in others only one is present.

I am unsure how to interpret these suffixes and what they represent in the context of the BANKSY embedding or feature space. I would be very grateful if you could provide some clarification on this. If my approach to finding differences among genes across neighbourhoods isn't appropriate, I would be happy to take suggestions.

Thank you very much for your time and assistance.

Regards,
Chiranjit

Code below:
test=readRDS("/cleaned/test.RDS")
test=JoinLayers(test)
test=NormalizeData(test)
test=FindVariableFeatures(test, nfeatures=4000)
hvg=VariableFeatures(test)
test=ScaleData(test, features=hvg)

group_name="sample"
dimx="CenterX_global_px"
dimy="CenterY_global_px"

test =RunBanksy(test, lambda = 0.8, assay = 'RNA', slot = 'data',
group = group_name, dimx = dimx, dimy = dimy,
k_geom = 15, split.scale = TRUE, use_agf=TRUE)

DefaultAssay(test)="BANKSY"

test = RunPCA(test, assay="BANKSY", features = rownames(test), npcs = 15)
test = RunHarmony(test, group.by.vars=group_name,
reduction.name='pca', reduction.save='harmony')

test = FindNeighbors(test, dims = 1:15, reduction = 'harmony')
test = FindClusters(test, resolution = 0.5, graph.name = 'BANKSY_snn')
test = FindClusters(test, resolution = 0.3, graph.name = 'BANKSY_snn')

markers=FindAllMarkers(test, assay="BANKSY", group.by="BANKSY_snn_res.0.3", logfc.threshold=0.25, min.pct = 0.2, only.pos=T)

[jleechung]

Hi @Chiranjit1504 sorry for the late reply --

As a first pass, I'd recommend running FindAllMarkers with assay set to your default assay (e.g. RNA) but with clusters determined by BANKSY clustering, something like this:

markers <- FindAllMarkers(test, assay = 'RNA', group.by="BANKSY_snn_res.0.3", logfc.threshold=0.25, min.pct = 0.2, only.pos=T)

This will find markers in your original feature set that are representative of each spatial domain BANKSY identifies, which makes it much easier to interpret.

The .m0 and .m1 suffixes represent the neighborhood-augmented feature sets that BANSKY computes. For instance:

• IL7R.m0 represents the mean expression of IL7R in the neighborhood of the cell -- this is relatively easy to interpret, meaning that for the cells in this cluster, the neighborhood average of IL7R is enriched vs. other clusters
• IL7R.m1 represents the spatial gradient of IL7R in the neighborhood of the index cell -- this is a higher order quantity, so less easy to interpret, but can roughly be interpreted as the spatial variation of that feature within the neighborhood

To your original point on finding differences between neighborhoods -- the BANKSY cluster labels already incorporate this neighborhood information, so switching back to the original gene set with assay = 'RNA' will find markers that are differentially expressed across neighborhoods.

Let me know if you have further questions,

[Chiranjit1504]

Thank you for the elaborate explanation @jleechung. Much appreciated.

Regards,
Chiranjit