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.
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:
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:
Finally, 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 4Parsing 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.000We 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.543nmi <- 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.675Visualise 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.363882 minssessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 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_2.15.0 AnnotationHub_3.15.0
#> [3] BiocFileCache_2.15.0 dbplyr_2.5.0
#> [5] spatialLIBD_1.18.0 cowplot_1.1.3
#> [7] scater_1.35.0 ggplot2_3.5.1
#> [9] harmony_1.2.3 Rcpp_1.0.13-1
#> [11] data.table_1.16.2 scran_1.35.0
#> [13] scuttle_1.17.0 Seurat_5.1.0
#> [15] SeuratObject_5.0.2 sp_2.1-4
#> [17] SpatialExperiment_1.17.0 SingleCellExperiment_1.29.1
#> [19] SummarizedExperiment_1.37.0 Biobase_2.67.0
#> [21] GenomicRanges_1.59.1 GenomeInfoDb_1.43.2
#> [23] IRanges_2.41.1 S4Vectors_0.45.2
#> [25] BiocGenerics_0.53.3 generics_0.1.3
#> [27] MatrixGenerics_1.19.0 matrixStats_1.4.1
#> [29] Banksy_1.3.0 BiocStyle_2.35.0
#>
#> loaded via a namespace (and not attached):
#> [1] bitops_1.0-9 spatstat.sparse_3.1-0 doParallel_1.0.17
#> [4] httr_1.4.7 RColorBrewer_1.1-3 tools_4.4.2
#> [7] sctransform_0.4.1 DT_0.33 utf8_1.2.4
#> [10] R6_2.5.1 lazyeval_0.2.2 uwot_0.2.2
#> [13] withr_3.0.2 gridExtra_2.3 progressr_0.15.1
#> [16] cli_3.6.3 spatstat.explore_3.3-3 fastDummies_1.7.4
#> [19] labeling_0.4.3 sass_0.4.9 spatstat.data_3.1-4
#> [22] ggridges_0.5.6 pbapply_1.7-2 Rsamtools_2.23.1
#> [25] dbscan_1.2-0 aricode_1.0.3 dichromat_2.0-0.1
#> [28] sessioninfo_1.2.2 parallelly_1.39.0 attempt_0.3.1
#> [31] maps_3.4.2.1 limma_3.63.2 pals_1.9
#> [34] RSQLite_2.3.8 BiocIO_1.17.1 ica_1.0-3
#> [37] spatstat.random_3.3-2 dplyr_1.1.4 Matrix_1.7-1
#> [40] ggbeeswarm_0.7.2 fansi_1.0.6 abind_1.4-8
#> [43] lifecycle_1.0.4 yaml_2.3.10 edgeR_4.5.0
#> [46] SparseArray_1.7.2 Rtsne_0.17 paletteer_1.6.0
#> [49] grid_4.4.2 blob_1.2.4 promises_1.3.2
#> [52] dqrng_0.4.1 crayon_1.5.3 miniUI_0.1.1.1
#> [55] lattice_0.22-6 beachmat_2.23.2 mapproj_1.2.11
#> [58] KEGGREST_1.47.0 magick_2.8.5 sys_3.4.3
#> [61] maketools_1.3.1 pillar_1.9.0 knitr_1.49
#> [64] metapod_1.15.0 rjson_0.2.23 future.apply_1.11.3
#> [67] codetools_0.2-20 leiden_0.4.3.1 glue_1.8.0
#> [70] spatstat.univar_3.1-1 vctrs_0.6.5 png_0.1-8
#> [73] spam_2.11-0 gtable_0.3.6 rematch2_2.1.2
#> [76] cachem_1.1.0 xfun_0.49 S4Arrays_1.7.1
#> [79] mime_0.12 survival_3.7-0 RcppHungarian_0.3
#> [82] iterators_1.0.14 fields_16.3 statmod_1.5.0
#> [85] bluster_1.17.0 fitdistrplus_1.2-1 ROCR_1.0-11
#> [88] nlme_3.1-166 bit64_4.5.2 filelock_1.0.3
#> [91] RcppAnnoy_0.0.22 bslib_0.8.0 irlba_2.3.5.1
#> [94] vipor_0.4.7 KernSmooth_2.23-24 colorspace_2.1-1
#> [97] DBI_1.2.3 tidyselect_1.2.1 bit_4.5.0
#> [100] compiler_4.4.2 curl_6.0.1 BiocNeighbors_2.1.1
#> [103] DelayedArray_0.33.2 plotly_4.10.4 rtracklayer_1.67.0
#> [106] scales_1.3.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-1 rmarkdown_2.29 benchmarkmeData_1.0.4
#> [115] RhpcBLASctl_0.23-42 XVector_0.47.0 htmltools_0.5.8.1
#> [118] pkgconfig_2.0.3 fastmap_1.2.0 rlang_1.1.4
#> [121] htmlwidgets_1.6.4 UCSC.utils_1.3.0 shiny_1.9.1
#> [124] farver_2.1.2 jquerylib_0.1.4 zoo_1.8-12
#> [127] jsonlite_1.8.9 BiocParallel_1.41.0 mclust_6.1.1
#> [130] config_0.3.2 RCurl_1.98-1.16 BiocSingular_1.23.0
#> [133] magrittr_2.0.3 GenomeInfoDbData_1.2.13 dotCall64_1.2
#> [136] patchwork_1.3.0 munsell_0.5.1 viridis_0.6.5
#> [139] reticulate_1.40.0 leidenAlg_1.1.4 stringi_1.8.4
#> [142] zlibbioc_1.52.0 MASS_7.3-61 plyr_1.8.9
#> [145] parallel_4.4.2 listenv_0.9.1 ggrepel_0.9.6
#> [148] deldir_2.0-4 Biostrings_2.75.1 sccore_1.0.5
#> [151] splines_4.4.2 tensor_1.5 locfit_1.5-9.10
#> [154] igraph_2.1.1 spatstat.geom_3.3-4 RcppHNSW_0.6.0
#> [157] buildtools_1.0.0 reshape2_1.4.4 ScaledMatrix_1.15.0
#> [160] XML_3.99-0.17 BiocVersion_3.21.1 evaluate_1.0.1
#> [163] golem_0.5.1 BiocManager_1.30.25 foreach_1.5.2
#> [166] httpuv_1.6.15 RANN_2.6.2 tidyr_1.3.1
#> [169] purrr_1.0.2 polyclip_1.10-7 benchmarkme_1.0.8
#> [172] future_1.34.0 scattermore_1.2 rsvd_1.0.5
#> [175] xtable_1.8-4 restfulr_0.0.15 RSpectra_0.16-2
#> [178] later_1.4.1 viridisLite_0.4.2 tibble_3.2.1
#> [181] GenomicAlignments_1.43.0 memoise_2.0.1 beeswarm_0.4.0
#> [184] AnnotationDbi_1.69.0 cluster_2.1.6 shinyWidgets_0.8.7
#> [187] globals_0.16.3