# Multi-sample analysis (10x Visium Human DLPFC)

Here, we demonstrate BANKSY analysis on 10x Visium data of the human dorsolateral prefrontal cortex from Maynard et al (2018). The data comprise 12 samples obtained from 3 subjects, with manual annotation of the layers in each sample. We will focus on 4 of the 12 samples from subject 3, demonstrating multi-sample analysis with BANKSY.

    library(Banksy)
    library(SummarizedExperiment)
    library(SpatialExperiment)
    library(Seurat)

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

# Loading the data

We fetch the data for all 12 DLPFC samples with the spatialLIBD package. This might take awhile.

    library(spatialLIBD)
    library(ExperimentHub)

    ehub <- ExperimentHub::ExperimentHub()
    spe <- spatialLIBD::fetch_data(type = "spe", eh = ehub)

After the download is completed, we trim the SpatialExperiment object, retaining only the counts and some metadata such as the sample identifier and pathology annotations. This saves some memory.

    imgData(spe) <- NULL
    assay(spe, "logcounts") <- NULL
    reducedDims(spe) <- NULL
    rowData(spe) <- NULL
    colData(spe) <- DataFrame(
        sample_id = spe$sample_id,
        clust_annotation = factor(
            addNA(spe$layer_guess_reordered_short),
            exclude = NULL, labels = seq(8)
        ),
        in_tissue = spe$in_tissue,
        row.names = colnames(spe)
    )
    invisible(gc())

Next, subset `spe` to samples from the last subject (samples `151673`, `151674`, `151675`, `151676`). This stores each sample in its own SpatialExperiment object, all placed in a list.

    sample_names <- as.character(151673:151676)
    spe_list <- lapply(sample_names, function(x) spe[, spe$sample_id == x])
    rm(spe)
    invisible(gc())

# Data preprocessing

Using Seurat, we perform basic normalisation of the data, and select the top 2000 highly variable features from each sample. Other methods for normalisation and feature selection may also be used. We take the union of these features for downstream analysis.

    #' Normalize data
    seu_list <- lapply(spe_list, function(x) {
        x <- as.Seurat(x, data = NULL)
        NormalizeData(x, scale.factor = 5000, normalization.method = 'RC')
    })

    #' Compute HVGs
    hvgs <- lapply(seu_list, function(x) {
        VariableFeatures(FindVariableFeatures(x, nfeatures = 2000))
    })
    hvgs <- Reduce(union, hvgs)

    #' Add data to SpatialExperiment and subset to HVGs
    aname <- "normcounts"
    spe_list <- Map(function(spe, seu) {
        assay(spe, aname) <- GetAssayData(seu)
        spe[hvgs,]
        }, spe_list, seu_list)
    rm(seu_list)
    invisible(gc())

# Running BANKSY

To run BANKSY across multiple samples, we first compute the BANKSY neighborhood feature matrices for each sample separately. We use `k_geom=6` corresponding to the first-order neighbors in 10x Visium assays (`k_geom=18` corresponding to first and second-order neighbors may also be used).

    compute_agf <- FALSE
    k_geom <- 6
    spe_list <- lapply(spe_list, computeBanksy, assay_name = aname,
                       compute_agf = compute_agf, k_geom = k_geom)

We then merge the samples to perform joint dimensional reduction and clustering:

    spe_joint <- do.call(cbind, spe_list)
    rm(spe_list)
    invisible(gc())

When running multi-sample BANKSY PCA, the `group` argument may be provided. This specifies the grouping variable for the cells or spots across the samples. Features belonging to cells or spots corresponding to each level of the grouping variable will be z-scaled separately. In this case, `sample_id` in `colData(spe_joint)` gives the grouping based on the sample of origin.

    lambda <- 0.2
    use_agf <- FALSE
    spe_joint <- runBanksyPCA(spe_joint, use_agf = use_agf, lambda = lambda, group = "sample_id", seed = 1000)

Run UMAP on the BANKSY embedding:

    spe_joint <- runBanksyUMAP(spe_joint, use_agf = use_agf, lambda = lambda, seed = 1000)

Finally, we obtain cluster labels for spots across all 4 samples. We use `connectClusters` for visual comparison of the manual annotations and BANKSY clusters.

    res <- 0.7
    spe_joint <- clusterBanksy(spe_joint, use_agf = use_agf, lambda = lambda, resolution = res, seed = 1000)
    cnm <- sprintf("clust_M%s_lam%s_k50_res%s", as.numeric(use_agf), lambda, res)
    spe_joint <- connectClusters(spe_joint)

Once joint clustering is performed, we split the samples into their own `SpatialExperiment` objects:

    spe_list <- lapply(sample_names, function(x) spe_joint[, spe_joint$sample_id == x])
    rm(spe_joint)
    invisible(gc())

As an optional step, we smooth the cluster labels of each sample separately. This can be done if smooth spatial domains are expected in the biological sample or tissue in question.

    spe_list <- lapply(spe_list, smoothLabels, cluster_names = cnm, k = 6L, verbose = FALSE)
    names(spe_list) <- paste0("sample_", sample_names)

The raw and smoothed cluster labels are stored in the `colData` slot of each `SingleCellExperiment` or `SpatialExperiment` object.

    #> DataFrame with 6 rows and 5 columns
    #>                      sample_id clust_annotation in_tissue
    #>                    <character>         <factor> <logical>
    #> AAACAAGTATCTCCCA-1      151673                3      TRUE
    #> AAACAATCTACTAGCA-1      151673                1      TRUE
    #> AAACACCAATAACTGC-1      151673                7      TRUE
    #> AAACAGAGCGACTCCT-1      151673                3      TRUE
    #> AAACAGCTTTCAGAAG-1      151673                5      TRUE
    #> AAACAGGGTCTATATT-1      151673                6      TRUE
    #>                    clust_M0_lam0.2_k50_res0.7 clust_M0_lam0.2_k50_res0.7_smooth
    #>                                      <factor>                          <factor>
    #> AAACAAGTATCTCCCA-1                          3                                 3
    #> AAACAATCTACTAGCA-1                          1                                 1
    #> AAACACCAATAACTGC-1                          7                                 7
    #> AAACAGAGCGACTCCT-1                          3                                 3
    #> AAACAGCTTTCAGAAG-1                          6                                 6
    #> AAACAGGGTCTATATT-1                          4                                 4

# Parsing BANKSY output

We can compare BANKSY clusters to pathology annotations using several cluster comparison measures such as the adjusted Rand index (ARI) or normalized mutual information (NMI) with `compareClusters`. The function computes the selected comparison measure between all pairs of cluster labels:

    compareClusters(spe_list$sample_151673, func = 'ARI')
    #>                                   clust_annotation clust_M0_lam0.2_k50_res0.7
    #> clust_annotation                             1.000                      0.548
    #> clust_M0_lam0.2_k50_res0.7                   0.548                      1.000
    #> clust_M0_lam0.2_k50_res0.7_smooth            0.566                      0.885
    #>                                   clust_M0_lam0.2_k50_res0.7_smooth
    #> clust_annotation                                              0.566
    #> clust_M0_lam0.2_k50_res0.7                                    0.885
    #> clust_M0_lam0.2_k50_res0.7_smooth                             1.000

We evaluate the ARI and NMI for each sample:

    ari <- sapply(spe_list, function(x) as.numeric(tail(compareClusters(x, func = "ARI")[, 1], n = 1)))
    ari
    #> sample_151673 sample_151674 sample_151675 sample_151676
    #>         0.566         0.575         0.529         0.543

    nmi <- sapply(spe_list, function(x) as.numeric(tail(compareClusters(x, func = "NMI")[, 1], n = 1)))
    nmi
    #> sample_151673 sample_151674 sample_151675 sample_151676
    #>         0.690         0.697         0.671         0.675

Visualise pathology annotation and BANKSY cluster on spatial coordinates with the ggspavis package:

    # Use scater:::.get_palette('tableau10medium')
    pal <- c(
        "#729ECE", "#FF9E4A", "#67BF5C", "#ED665D", "#AD8BC9",
        "#A8786E", "#ED97CA", "#A2A2A2", "#CDCC5D", "#6DCCDA"
    )

    plot_bank <- lapply(spe_list, function(x) {
        df <- cbind.data.frame(
            clust=colData(x)[[sprintf("%s_smooth", cnm)]], spatialCoords(x))
        ggplot(df, aes(x=pxl_row_in_fullres, y=pxl_col_in_fullres, col=clust)) +
            geom_point(size = 0.5) +
            scale_color_manual(values = pal) +
            theme_classic() +
            theme(
                legend.position = "none",
                axis.text.x=element_blank(),
                axis.text.y=element_blank(),
                axis.ticks=element_blank(),
                axis.title.x=element_blank(),
                axis.title.y=element_blank()) +
            labs(title = "BANKSY clusters") +
            coord_equal()
    })

    plot_anno <- lapply(spe_list, function(x) {
        df <- cbind.data.frame(
            clust=colData(x)[['clust_annotation']], spatialCoords(x))
        ggplot(df, aes(x=pxl_row_in_fullres, y=pxl_col_in_fullres, col=clust)) +
            geom_point(size = 0.5) +
            scale_color_manual(values = pal) +
            theme_classic() +
            theme(
                legend.position = "none",
                axis.text.x=element_blank(),
                axis.text.y=element_blank(),
                axis.ticks=element_blank(),
                axis.title.x=element_blank(),
                axis.title.y=element_blank()) +
            labs(title = sprintf("Sample %s", x$sample_id[1])) +
            coord_equal()
    })

    plot_list <- c(plot_anno, plot_bank)

    plot_grid(plotlist = plot_list, ncol = 4, byrow = TRUE)

Visualize joint UMAPs for each sample:

    umap_bank <- lapply(spe_list, function(x) {
        plotReducedDim(x,
            "UMAP_M0_lam0.2",
            colour_by = sprintf("%s_smooth", cnm),
            point_size = 0.5
        ) +
            theme(legend.position = "none") +
            labs(title = "BANKSY clusters")
    })

    umap_anno <- lapply(spe_list, function(x) {
        plotReducedDim(x,
            "UMAP_M0_lam0.2",
            colour_by = "clust_annotation",
            point_size = 0.5
        ) +
            theme(legend.position = "none") +
            labs(title = sprintf("Sample %s", x$sample_id[1]))
    })

    umap_list <- c(umap_anno, umap_bank)

    plot_grid(plotlist = umap_list, ncol = 4, byrow = TRUE)

# Session information

Vignette runtime:

    #> Time difference of 1.402369 mins

    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