# Spatial data integration with Harmony (10x Visium Human DLPFC)

#### 4 June 2025

# Contents

Here, we demonstrate how BANKSY can be used with Harmony for integrating multiple spatial omics datasets in the presence of strong batch effects. We use 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.

    library(Banksy)
    library(SummarizedExperiment)
    library(SpatialExperiment)
    library(Seurat)
    library(scran)
    library(data.table)
    library(harmony)

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

    SEED <- 1000

# 1 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.

    #' Remove NA spots
    na_id <- which(is.na(spe$layer_guess_reordered_short))
    spe <- spe[, -na_id]

    #' Trim
    imgData(spe) <- NULL
    assay(spe, "logcounts") <- NULL
    reducedDims(spe) <- NULL
    rowData(spe) <- NULL
    colData(spe) <- DataFrame(
        sample_id = spe$sample_id,
        subject_id = factor(spe$sample_id, labels = rep(paste0("Subject", 1:3), each = 4)),
        clust_annotation = factor(as.numeric(spe$layer_guess_reordered_short)),
        in_tissue = spe$in_tissue,
        row.names = colnames(spe)
    )
    colnames(spe) <- paste0(colnames(spe), "_", spe$sample_id)
    invisible(gc())

We analyse the first sample of each subject due to vignette runtime constraints.

    spe <- spe[, spe$sample_id %in% c("151507", "151669", "151673")]
    sample_names <- unique(spe$sample_id)

Next, stagger the spatial coordinates across the samples so that spots from different samples do not overlap.

    #' Stagger spatial coordinates
    locs <- spatialCoords(spe)
    locs <- cbind(locs, sample_id = factor(spe$sample_id))
    locs_dt <- data.table(locs)
    colnames(locs_dt) <- c("sdimx", "sdimy", "group")
    locs_dt[, sdimx := sdimx - min(sdimx), by = group]
    global_max <- max(locs_dt$sdimx) * 1.5
    locs_dt[, sdimx := sdimx + group * global_max]
    locs <- as.matrix(locs_dt[, 1:2])
    rownames(locs) <- colnames(spe)
    spatialCoords(spe) <- locs

# 2 Data preprocessing

Find highly variable features and normalize counts. Here we use Seurat, but other methods may also be used (e.g. `scran::getTopHVGs`).

    #' Get HVGs
    seu <- as.Seurat(spe, data = NULL)
    seu <- FindVariableFeatures(seu, nfeatures = 2000)

    #' Normalize data
    scale_factor <- median(colSums(assay(spe, "counts")))
    seu <- NormalizeData(seu, scale.factor = scale_factor, normalization.method = "RC")

    #' Add data to SpatialExperiment and subset to HVGs
    aname <- "normcounts"
    assay(spe, aname) <- GetAssayData(seu)
    spe <- spe[VariableFeatures(seu),]

# 3 Running BANKSY

Compute BANKSY neighborhood matrices. We use `k_geom=18` corresponding to first and second-order neighbors in 10x Visium.

    compute_agf <- TRUE
    k_geom <- 18
    spe <- computeBanksy(spe, assay_name = aname, compute_agf = compute_agf, k_geom = k_geom)

Run PCA on the BANKSY matrix:

    lambda <- 0.2
    npcs <- 20
    use_agf <- TRUE
    spe <- runBanksyPCA(spe, use_agf = use_agf, lambda = lambda, npcs = npcs, seed = SEED)

# 4 Run Harmony on BANKSY’s embedding

We run Harmony on the PCs of the BANKSY matrix:

    set.seed(SEED)
    harmony_embedding <- HarmonyMatrix(
        data_mat = reducedDim(spe, "PCA_M1_lam0.2"),
        meta_data = colData(spe),
        vars_use = c("sample_id", "subject_id"),
        do_pca = FALSE,
        max.iter.harmony = 20,
        verbose = FALSE
    )
    reducedDim(spe, "Harmony_BANKSY") <- harmony_embedding

Next, run UMAP on the ‘raw’ and Harmony corrected PCA embeddings:

    spe <- runBanksyUMAP(spe, use_agf = TRUE, lambda = lambda, npcs = npcs)
    spe <- runBanksyUMAP(spe, dimred = "Harmony_BANKSY")

Visualize the UMAPs annotated by subject ID:

    plot_grid(
        plotReducedDim(spe, "UMAP_M1_lam0.2",
                       point_size = 0.1,
                       point_alpha = 0.5,
                       color_by = "subject_id") +
            theme(legend.position = "none"),
        plotReducedDim(spe, "UMAP_Harmony_BANKSY",
                       point_size = 0.1,
                       point_alpha = 0.5,
                       color_by = "subject_id") +
            theme(legend.title = element_blank()) +
            guides(colour = guide_legend(override.aes = list(size = 5, alpha = 1))),
        nrow = 1,
        rel_widths = c(1, 1.2)
    )

Cluster the Harmony corrected PCA embedding:

    spe <- clusterBanksy(spe, dimred = "Harmony_BANKSY", resolution = 0.55, seed = SEED)
    spe <- connectClusters(spe, map_to = "clust_annotation")

Generate spatial plots:

    cnm <- clusterNames(spe)[2]
    spatial_plots <- lapply(sample_names, function(snm) {
        x <- spe[, spe$sample_id == snm]
        ari <- aricode::ARI(x$clust_annotation, colData(x)[, cnm])

        df <- cbind.data.frame(clust=colData(x)[[cnm]], spatialCoords(x))
        ggplot(df, aes(x=sdimy, y=sdimx, col=clust)) +
            geom_point(size = 0.5) +
            scale_color_manual(values = pals::kelly()[-1]) +
            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 - ARI: %s", snm, round(ari, 3))) +
            coord_equal()
    })

    plot_grid(plotlist = spatial_plots, ncol = 3, byrow = FALSE)

# 5 Session information

Vignette runtime:

    #> Time difference of 4.152628 mins

    sessionInfo()
    #> R version 4.5.0 (2025-04-11 ucrt)
    #> Platform: x86_64-w64-mingw32/x64
    #> Running under: Windows Server 2022 x64 (build 20348)
    #>
    #> Matrix products: default
    #>   LAPACK version 3.12.1
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.utf8
    #> [3] LC_MONETARY=English_United States.utf8
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.utf8
    #>
    #> time zone: America/New_York
    #> tzcode source: internal
    #>
    #> attached base packages:
    #> [1] stats4    stats     graphics  grDevices utils     datasets  methods
    #> [8] base
    #>
    #> other attached packages:
    #>  [1] ExperimentHub_2.99.5        AnnotationHub_3.99.5
    #>  [3] BiocFileCache_2.99.5        dbplyr_2.5.0
    #>  [5] spatialLIBD_1.21.5          cowplot_1.1.3
    #>  [7] scater_1.37.0               ggplot2_3.5.2
    #>  [9] harmony_1.2.3               Rcpp_1.0.14
    #> [11] data.table_1.17.4           scran_1.37.0
    #> [13] scuttle_1.19.0              Seurat_5.3.0
    #> [15] SeuratObject_5.1.0          sp_2.2-0
    #> [17] SpatialExperiment_1.19.1    SingleCellExperiment_1.31.0
    #> [19] SummarizedExperiment_1.39.0 Biobase_2.69.0
    #> [21] GenomicRanges_1.61.0        GenomeInfoDb_1.45.4
    #> [23] IRanges_2.43.0              S4Vectors_0.47.0
    #> [25] BiocGenerics_0.55.0         generics_0.1.4
    #> [27] MatrixGenerics_1.21.0       matrixStats_1.5.0
    #> [29] Banksy_1.5.0                BiocStyle_2.37.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.0
    #>   [7] sctransform_0.4.2        DT_0.33                  R6_2.6.1
    #>  [10] lazyeval_0.2.2           uwot_0.2.3               GetoptLong_1.0.5
    #>  [13] withr_3.0.2              gridExtra_2.3            progressr_0.15.1
    #>  [16] cli_3.6.5                spatstat.explore_3.4-3   fastDummies_1.7.5
    #>  [19] labeling_0.4.3           sass_0.4.10              spatstat.data_3.1-6
    #>  [22] ggridges_0.5.6           pbapply_1.7-2            Rsamtools_2.25.0
    #>  [25] dbscan_1.2.2             aricode_1.0.3            dichromat_2.0-0.1
    #>  [28] sessioninfo_1.2.3        parallelly_1.45.0        attempt_0.3.1
    #>  [31] maps_3.4.3               limma_3.65.1             pals_1.10
    #>  [34] RSQLite_2.4.0            BiocIO_1.19.0            shape_1.4.6.1
    #>  [37] ica_1.0-3                spatstat.random_3.4-1    dplyr_1.1.4
    #>  [40] Matrix_1.7-3             ggbeeswarm_0.7.2         abind_1.4-8
    #>  [43] lifecycle_1.0.4          yaml_2.3.10              edgeR_4.7.2
    #>  [46] SparseArray_1.9.0        Rtsne_0.17               paletteer_1.6.0
    #>  [49] grid_4.5.0               blob_1.2.4               promises_1.3.3
    #>  [52] dqrng_0.4.1              crayon_1.5.3             miniUI_0.1.2
    #>  [55] lattice_0.22-7           beachmat_2.25.1          mapproj_1.2.12
    #>  [58] KEGGREST_1.49.0          magick_2.8.6             pillar_1.10.2
    #>  [61] knitr_1.50               ComplexHeatmap_2.25.0    metapod_1.17.0
    #>  [64] rjson_0.2.23             future.apply_1.11.3      codetools_0.2-20
    #>  [67] glue_1.8.0               spatstat.univar_3.1-3    vctrs_0.6.5
    #>  [70] png_0.1-8                spam_2.11-1              gtable_0.3.6
    #>  [73] rematch2_2.1.2           cachem_1.1.0             xfun_0.52
    #>  [76] S4Arrays_1.9.1           mime_0.13                survival_3.8-3
    #>  [79] RcppHungarian_0.3        iterators_1.0.14         tinytex_0.57
    #>  [82] statmod_1.5.0            bluster_1.19.0           fitdistrplus_1.2-2
    #>  [85] ROCR_1.0-11              nlme_3.1-168             bit64_4.6.0-1
    #>  [88] filelock_1.0.3           RcppAnnoy_0.0.22         bslib_0.9.0
    #>  [91] irlba_2.3.5.1            vipor_0.4.7              KernSmooth_2.23-26
    #>  [94] colorspace_2.1-1         DBI_1.2.3                tidyselect_1.2.1
    #>  [97] bit_4.6.0                compiler_4.5.0           curl_6.2.3
    #> [100] httr2_1.1.2              BiocNeighbors_2.3.1      DelayedArray_0.35.1
    #> [103] plotly_4.10.4            rtracklayer_1.69.0       bookdown_0.43
    #> [106] scales_1.4.0             lmtest_0.9-40            rappdirs_0.3.3
    #> [109] stringr_1.5.1            digest_0.6.37            goftest_1.2-3
    #> [112] spatstat.utils_3.1-4     rmarkdown_2.29           benchmarkmeData_1.0.4
    #> [115] RhpcBLASctl_0.23-42      XVector_0.49.0           htmltools_0.5.8.1
    #> [118] pkgconfig_2.0.3          fastmap_1.2.0            GlobalOptions_0.1.2
    #> [121] rlang_1.1.6              htmlwidgets_1.6.4        UCSC.utils_1.5.0
    #> [124] shiny_1.10.0             farver_2.1.2             jquerylib_0.1.4
    #> [127] zoo_1.8-14               jsonlite_2.0.0           BiocParallel_1.43.3
    #> [130] mclust_6.1.1             config_0.3.2             RCurl_1.98-1.17
    #> [133] BiocSingular_1.25.0      magrittr_2.0.3           dotCall64_1.2
    #> [136] patchwork_1.3.0          viridis_0.6.5            reticulate_1.42.0
    #> [139] leidenAlg_1.1.5          stringi_1.8.7            MASS_7.3-65
    #> [142] plyr_1.8.9               parallel_4.5.0           listenv_0.9.1
    #> [145] ggrepel_0.9.6            deldir_2.0-4             Biostrings_2.77.1
    #> [148] sccore_1.0.6             splines_4.5.0            tensor_1.5
    #> [151] circlize_0.4.16          locfit_1.5-9.12          igraph_2.1.4
    #> [154] spatstat.geom_3.4-1      RcppHNSW_0.6.0           reshape2_1.4.4
    #> [157] ScaledMatrix_1.17.0      XML_3.99-0.18            BiocVersion_3.22.0
    #> [160] evaluate_1.0.3           golem_0.5.1              BiocManager_1.30.25
    #> [163] foreach_1.5.2            httpuv_1.6.16            RANN_2.6.2
    #> [166] tidyr_1.3.1              purrr_1.0.4              polyclip_1.10-7
    #> [169] benchmarkme_1.0.8        clue_0.3-66              future_1.49.0
    #> [172] scattermore_1.2          rsvd_1.0.5               xtable_1.8-4
    #> [175] restfulr_0.0.15          RSpectra_0.16-2          later_1.4.2
    #> [178] viridisLite_0.4.2        tibble_3.2.1             GenomicAlignments_1.45.0
    #> [181] memoise_2.0.1            beeswarm_0.4.0           AnnotationDbi_1.71.0
    #> [184] cluster_2.1.8.1          shinyWidgets_0.9.0       globals_0.18.0