# Parameter selection (VeraFISH Mouse Hippocampus)

Here, we demonstrate a grid search of clustering parameters with a mouse hippocampus VeraFISH dataset. BANKSY currently provides four algorithms for clustering the BANKSY matrix with clusterBanksy: Leiden (default), Louvain, k-means, and model-based clustering. In this vignette, we run only Leiden clustering. See `?clusterBanksy` for more details on the parameters for different clustering methods.

# Loading the data

The dataset comprises gene expression for 10,944 cells and 120 genes in 2 spatial dimensions. See `?Banksy::hippocampus` for more details.

    # Load libs
    library(Banksy)

    library(SummarizedExperiment)
    library(SpatialExperiment)
    library(scuttle)

    library(scater)
    library(cowplot)
    library(ggplot2)

    # Load data
    data(hippocampus)
    gcm <- hippocampus$expression
    locs <- as.matrix(hippocampus$locations)

Here, `gcm` is a gene by cell matrix, and `locs` is a matrix specifying the coordinates of the centroid for each cell.

    head(gcm[,1:5])
    #>         cell_1276 cell_8890 cell_691 cell_396 cell_9818
    #> Sparcl1        45         0       11       22         0
    #> Slc1a2         17         0        6        5         0
    #> Map            10         0       12       16         0
    #> Sqstm1         26         0        0        2         0
    #> Atp1a2          0         0        4        3         0
    #> Tnc             0         0        0        0         0
    head(locs)
    #>                 sdimx    sdimy
    #> cell_1276  -13372.899 15776.37
    #> cell_8890    8941.101 15866.37
    #> cell_691   -14882.899 15896.37
    #> cell_396   -15492.899 15835.37
    #> cell_9818   11308.101 15846.37
    #> cell_11310  14894.101 15810.37

Initialize a SpatialExperiment object and perform basic quality control. We keep cells with total transcript count within the 5th and 98th percentile:

    se <- SpatialExperiment(assay = list(counts = gcm), spatialCoords = locs)
    colData(se) <- cbind(colData(se), spatialCoords(se))

    # QC based on total counts
    qcstats <- perCellQCMetrics(se)
    thres <- quantile(qcstats$total, c(0.05, 0.98))
    keep <- (qcstats$total > thres[1]) & (qcstats$total < thres[2])
    se <- se[, keep]

Next, perform normalization of the data.

    # Normalization to mean library size
    se <- computeLibraryFactors(se)
    aname <- "normcounts"
    assay(se, aname) <- normalizeCounts(se, log = FALSE)

# Parameters

BANKSY has a few key parameters. We describe these below.

## AGF usage

For characterising neighborhoods, BANKSY computes the weighted neighborhood mean (`H_0`) and the azimuthal Gabor filter (`H_1`), which estimates gene expression gradients. Setting `compute_agf=TRUE` computes both `H_0` and `H_1`.

## k-geometric

`k_geom` specifies the number of neighbors used to compute each `H_m` for `m=0,1`. If a single value is specified, the same `k_geom` will be used for each feature matrix. Alternatively, multiple values of `k_geom` can be provided for each feature matrix. Here, we use `k_geom[1]=15` and `k_geom[2]=30` for `H_0` and `H_1` respectively. More neighbors are used to compute gradients.

> For datasets generated using Visium v1/v2, use `k_geom=18` (or `k_geom <- c(18, 18)` if `compute_agf = TRUE`), since that corresponds to taking as neighbourhood two concentric rings of spots around each spot.

We compute the neighborhood feature matrices using normalized expression (`normcounts` in the `se` object).

    k_geom <- c(15, 30)
    se <- computeBanksy(se, assay_name = aname, compute_agf = TRUE, k_geom = k_geom)
    #> Computing neighbors...
    #> Spatial mode is kNN_median
    #> Parameters: k_geom=15
    #> Done
    #> Computing neighbors...
    #> Spatial mode is kNN_median
    #> Parameters: k_geom=30
    #> Done
    #> Done
    #> Centering
    #> Done

`computeBanksy` populates the `assays` slot with `H_0` and `H_1` in this instance:

    se
    #> class: SpatialExperiment
    #> dim: 120 10205
    #> metadata(1): BANKSY_params
    #> assays(4): counts normcounts H0 H1
    #> rownames(120): Sparcl1 Slc1a2 ... Notch3 Egfr
    #> rowData names(0):
    #> colnames(10205): cell_1276 cell_691 ... cell_11635 cell_10849
    #> colData names(4): sample_id sdimx sdimy sizeFactor
    #> reducedDimNames(0):
    #> mainExpName: NULL
    #> altExpNames(0):
    #> spatialCoords names(2) : sdimx sdimy
    #> imgData names(1): sample_id

## lambda

The `lambda` parameter is a mixing parameter in `[0,1]` which determines how much spatial information is incorporated for downstream analysis. With smaller values of `lambda`, BANKY operates in cell-typing mode, while at higher levels of `lambda`, BANKSY operates in domain-finding mode. As a starting point, we recommend `lambda=0.2` for cell-typing and `lambda=0.8` for zone-finding, except for datasets generated using the Visium v1/v2 technology, for which we recommend `lambda=0.2` for domain finding. See the note in the tutorial on the main page for more info.

Here, we run `lambda=0` which corresponds to non-spatial clustering, and `lambda=0.2` for spatially-informed cell-typing. We compute PCs with and without the AGF (`H_1`).

    lambda <- c(0, 0.2)
    se <- runBanksyPCA(se, use_agf = c(FALSE, TRUE), lambda = lambda, seed = 1000)
    #> Using seed=1000
    #> Using seed=1000
    #> Using seed=1000
    #> Using seed=1000

`runBanksyPCA` populates the `reducedDims` slot, with each combination of `use_agf` and `lambda` provided.

    reducedDimNames(se)
    #> [1] "PCA_M0_lam0"   "PCA_M0_lam0.2" "PCA_M1_lam0"   "PCA_M1_lam0.2"

## Clustering parameters

Next, we cluster the BANKSY embedding with Leiden graph-based clustering. This admits two parameters: `k_neighbors` and `resolution`. `k_neighbors` determines the number of k nearest neighbors used to construct the shared nearest neighbors graph. Leiden clustering is then performed on the resultant graph with resolution `resolution`. For reproducibiltiy we set a seed for each parameter combination.

    k <- 50
    res <- 1
    se <- clusterBanksy(se, use_agf = c(FALSE, TRUE), lambda = lambda, k_neighbors = k, resolution = res, seed = 1000)
    #> Using seed=1000
    #> Using seed=1000
    #> Using seed=1000
    #> Using seed=1000

`clusterBanksy` populates `colData(se)` with cluster labels:

    colnames(colData(se))
    #> [1] "sample_id"                "sdimx"
    #> [3] "sdimy"                    "sizeFactor"
    #> [5] "clust_M0_lam0_k50_res1"   "clust_M0_lam0.2_k50_res1"
    #> [7] "clust_M1_lam0_k50_res1"   "clust_M1_lam0.2_k50_res1"

# Comparing cluster results

To compare clustering runs visually, different runs can be relabeled to minimise their differences with `connectClusters`:

    se <- connectClusters(se)
    #> clust_M1_lam0_k50_res1 --> clust_M0_lam0_k50_res1
    #> clust_M0_lam0.2_k50_res1 --> clust_M1_lam0_k50_res1
    #> clust_M1_lam0.2_k50_res1 --> clust_M0_lam0.2_k50_res1

Visualise spatial coordinates with cluster labels.

    cnames <- colnames(colData(se))
    cnames <- cnames[grep("^clust", cnames)]
    cplots <- lapply(cnames, function(cnm) {
        plotColData(se, x = "sdimx", y = "sdimy", point_size = 0.1, colour_by = cnm) +
            coord_equal() +
            labs(title = cnm) +
            theme(legend.title = element_blank()) +
            guides(colour = guide_legend(override.aes = list(size = 2)))
    })

    plot_grid(plotlist = cplots, ncol = 2)

Compare all cluster outputs with `compareClusters`. This function computes pairwise cluster comparison metrics between the clusters in `colData(se)` based on adjusted Rand index (ARI):

    compareClusters(se, func = "ARI")
    #>                          clust_M0_lam0_k50_res1 clust_M0_lam0.2_k50_res1
    #> clust_M0_lam0_k50_res1                    1.000                     0.67
    #> clust_M0_lam0.2_k50_res1                  0.670                     1.00
    #> clust_M1_lam0_k50_res1                    1.000                     0.67
    #> clust_M1_lam0.2_k50_res1                  0.747                     0.87
    #>                          clust_M1_lam0_k50_res1 clust_M1_lam0.2_k50_res1
    #> clust_M0_lam0_k50_res1                    1.000                    0.747
    #> clust_M0_lam0.2_k50_res1                  0.670                    0.870
    #> clust_M1_lam0_k50_res1                    1.000                    0.747
    #> clust_M1_lam0.2_k50_res1                  0.747                    1.000

or normalized mutual information (NMI):

    compareClusters(se, func = "NMI")
    #>                          clust_M0_lam0_k50_res1 clust_M0_lam0.2_k50_res1
    #> clust_M0_lam0_k50_res1                    1.000                    0.741
    #> clust_M0_lam0.2_k50_res1                  0.741                    1.000
    #> clust_M1_lam0_k50_res1                    1.000                    0.741
    #> clust_M1_lam0.2_k50_res1                  0.782                    0.915
    #>                          clust_M1_lam0_k50_res1 clust_M1_lam0.2_k50_res1
    #> clust_M0_lam0_k50_res1                    1.000                    0.782
    #> clust_M0_lam0.2_k50_res1                  0.741                    0.915
    #> clust_M1_lam0_k50_res1                    1.000                    0.782
    #> clust_M1_lam0.2_k50_res1                  0.782                    1.000

See `?compareClusters` for the full list of comparison measures.

# Session information

Vignette runtime:

    #> Time difference of 33.14853 secs

    sessionInfo()
    #> R version 4.5.2 (2025-10-31)
    #> Platform: x86_64-pc-linux-gnu
    #> Running under: Ubuntu 24.04.3 LTS
    #>
    #> Matrix products: default
    #> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
    #> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
    #>
    #> locale:
    #>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
    #>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C
    #>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
    #>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
    #>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    #> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
    #>
    #> time zone: Etc/UTC
    #> tzcode source: system (glibc)
    #>
    #> attached base packages:
    #> [1] stats4    stats     graphics  grDevices utils     datasets  methods
    #> [8] base
    #>
    #> other attached packages:
    #>  [1] ExperimentHub_3.1.0         AnnotationHub_4.1.0
    #>  [3] BiocFileCache_3.1.0         dbplyr_2.5.1
    #>  [5] spatialLIBD_1.23.0          cowplot_1.2.0
    #>  [7] scater_1.39.1               ggplot2_4.0.1
    #>  [9] harmony_1.2.4               Rcpp_1.1.0.8.2
    #> [11] data.table_1.18.0           scran_1.39.0
    #> [13] scuttle_1.21.0              Seurat_5.4.0
    #> [15] SeuratObject_5.3.0          sp_2.2-0
    #> [17] SpatialExperiment_1.21.0    SingleCellExperiment_1.33.0
    #> [19] SummarizedExperiment_1.41.0 Biobase_2.71.0
    #> [21] GenomicRanges_1.63.1        Seqinfo_1.1.0
    #> [23] IRanges_2.45.0              S4Vectors_0.49.0
    #> [25] BiocGenerics_0.57.0         generics_0.1.4
    #> [27] MatrixGenerics_1.23.0       matrixStats_1.5.0
    #> [29] Banksy_1.7.0                BiocStyle_2.39.0
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] bitops_1.0-9             spatstat.sparse_3.1-0    httr_1.4.7
    #>   [4] RColorBrewer_1.1-3       doParallel_1.0.17        tools_4.5.2
    #>   [7] sctransform_0.4.2        DT_0.34.0                R6_2.6.1
    #>  [10] lazyeval_0.2.2           uwot_0.2.4               GetoptLong_1.1.0
    #>  [13] withr_3.0.2              gridExtra_2.3            progressr_0.18.0
    #>  [16] cli_3.6.5                spatstat.explore_3.6-0   fastDummies_1.7.5
    #>  [19] labeling_0.4.3           sass_0.4.10              S7_0.2.1
    #>  [22] spatstat.data_3.1-9      ggridges_0.5.7           pbapply_1.7-4
    #>  [25] Rsamtools_2.27.0         dbscan_1.2.4             aricode_1.0.3
    #>  [28] dichromat_2.0-0.1        sessioninfo_1.2.3        parallelly_1.46.1
    #>  [31] attempt_0.3.1            maps_3.4.3               limma_3.67.0
    #>  [34] pals_1.10                RSQLite_2.4.5            BiocIO_1.21.0
    #>  [37] shape_1.4.6.1            ica_1.0-3                spatstat.random_3.4-3
    #>  [40] dplyr_1.1.4              Matrix_1.7-4             ggbeeswarm_0.7.3
    #>  [43] abind_1.4-8              lifecycle_1.0.5          yaml_2.3.12
    #>  [46] edgeR_4.9.2              SparseArray_1.11.10      Rtsne_0.17
    #>  [49] paletteer_1.7.0          grid_4.5.2               blob_1.2.4
    #>  [52] promises_1.5.0           dqrng_0.4.1              crayon_1.5.3
    #>  [55] miniUI_0.1.2             lattice_0.22-7           beachmat_2.27.1
    #>  [58] cigarillo_1.1.0          mapproj_1.2.12           KEGGREST_1.51.1
    #>  [61] magick_2.9.0             sys_3.4.3                maketools_1.3.2
    #>  [64] pillar_1.11.1            knitr_1.51               ComplexHeatmap_2.27.0
    #>  [67] metapod_1.19.1           rjson_0.2.23             future.apply_1.20.1
    #>  [70] codetools_0.2-20         glue_1.8.0               spatstat.univar_3.1-5
    #>  [73] vctrs_0.6.5              png_0.1-8                spam_2.11-3
    #>  [76] gtable_0.3.6             rematch2_2.1.2           cachem_1.1.0
    #>  [79] xfun_0.55                S4Arrays_1.11.1          mime_0.13
    #>  [82] survival_3.8-3           RcppHungarian_0.3        iterators_1.0.14
    #>  [85] statmod_1.5.1            bluster_1.21.0           fitdistrplus_1.2-4
    #>  [88] ROCR_1.0-11              nlme_3.1-168             bit64_4.6.0-1
    #>  [91] filelock_1.0.3           RcppAnnoy_0.0.22         bslib_0.9.0
    #>  [94] irlba_2.3.5.1            vipor_0.4.7              KernSmooth_2.23-26
    #>  [97] otel_0.2.0               colorspace_2.1-2         DBI_1.2.3
    #> [100] tidyselect_1.2.1         bit_4.6.0                compiler_4.5.2
    #> [103] curl_7.0.0               httr2_1.2.2              BiocNeighbors_2.5.0
    #> [106] DelayedArray_0.37.0      plotly_4.11.0            rtracklayer_1.71.3
    #> [109] scales_1.4.0             lmtest_0.9-40            rappdirs_0.3.3
    #> [112] stringr_1.6.0            digest_0.6.39            goftest_1.2-3
    #> [115] spatstat.utils_3.2-0     rmarkdown_2.30           benchmarkmeData_1.0.4
    #> [118] RhpcBLASctl_0.23-42      XVector_0.51.0           htmltools_0.5.9
    #> [121] pkgconfig_2.0.3          fastmap_1.2.0            GlobalOptions_0.1.3
    #> [124] rlang_1.1.7              htmlwidgets_1.6.4        shiny_1.12.1
    #> [127] farver_2.1.2             jquerylib_0.1.4          zoo_1.8-15
    #> [130] jsonlite_2.0.0           BiocParallel_1.45.0      mclust_6.1.2
    #> [133] config_0.3.2             RCurl_1.98-1.17          BiocSingular_1.27.1
    #> [136] magrittr_2.0.4           dotCall64_1.2            patchwork_1.3.2
    #> [139] viridis_0.6.5            reticulate_1.44.1        leidenAlg_1.1.5
    #> [142] stringi_1.8.7            MASS_7.3-65              plyr_1.8.9
    #> [145] parallel_4.5.2           listenv_0.10.0           ggrepel_0.9.6
    #> [148] deldir_2.0-4             Biostrings_2.79.4        sccore_1.0.6
    #> [151] splines_4.5.2            tensor_1.5.1             circlize_0.4.17
    #> [154] locfit_1.5-9.12          igraph_2.2.1             spatstat.geom_3.6-1
    #> [157] RcppHNSW_0.6.0           buildtools_1.0.0         reshape2_1.4.5
    #> [160] ScaledMatrix_1.19.0      XML_3.99-0.20            BiocVersion_3.23.1
    #> [163] evaluate_1.0.5           golem_0.5.1              BiocManager_1.30.27
    #> [166] foreach_1.5.2            httpuv_1.6.16            RANN_2.6.2
    #> [169] tidyr_1.3.2              purrr_1.2.1              polyclip_1.10-7
    #> [172] benchmarkme_1.0.8        clue_0.3-66              future_1.68.0
    #> [175] scattermore_1.2          rsvd_1.0.5               xtable_1.8-4
    #> [178] restfulr_0.0.16          RSpectra_0.16-2          later_1.4.5
    #> [181] viridisLite_0.4.2        tibble_3.3.0             GenomicAlignments_1.47.0
    #> [184] memoise_2.0.1            beeswarm_0.4.0           AnnotationDbi_1.73.0
    #> [187] cluster_2.1.8.1          shinyWidgets_0.9.0       globals_0.18.0