## ----global.options, include = FALSE------------------------------------------ knitr::opts_knit$set( collapse = TRUE, comment = "#>", fig.align = 'center' ) knitr::opts_chunk$set( out.extra = 'style="display:block; margin:auto;"', warning = FALSE, message = FALSE, # Compact, HTML-friendly rendering. JPEG (rather than PNG) is the # right choice for the H&E histology overlays — photographic raster # compresses ~5x better as JPEG at quality 80, keeping the rendered # vignette under Bioconductor's 5 MB per-file size guideline. dev = "jpeg", dev.args = list(quality = 80), dpi = 72, fig.width = 5, fig.height = 4, fig.retina = 1 ) ## ----eval = FALSE------------------------------------------------------------- # if (!require("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("SpaceMarkers") ## ----message=FALSE, warning=FALSE--------------------------------------------- library(SpaceMarkers) ## ----------------------------------------------------------------------------- data_dir <- "visiumBrCA" unlink(file.path(data_dir), recursive = TRUE) dir.create(data_dir, showWarnings = FALSE) main_10xlink <- "https://cf.10xgenomics.com/samples/spatial-exp/1.3.0" counts_folder <- "Visium_Human_Breast_Cancer" counts_file <- "Visium_Human_Breast_Cancer_filtered_feature_bc_matrix.h5" counts_url <- paste(c(main_10xlink, counts_folder, counts_file), collapse = "/") sp_folder <- "Visium_Human_Breast_Cancer" sp_file <- "Visium_Human_Breast_Cancer_spatial.tar.gz" sp_url <- paste(c(main_10xlink, sp_folder, sp_file), collapse = "/") ## ----------------------------------------------------------------------------- bfc <- BiocFileCache::BiocFileCache(ask = FALSE) counts_cache <- BiocFileCache::bfcrpath(bfc, counts_url) sp_cache <- BiocFileCache::bfcrpath(bfc, sp_url) file.copy(counts_cache, file.path(data_dir, counts_file), overwrite = TRUE) untar(sp_cache, exdir = data_dir) ## ----message=FALSE, warning=FALSE--------------------------------------------- cogaps_result <- readRDS(system.file("extdata", "CoGAPS_result.rds", package = "SpaceMarkers", mustWork = TRUE)) sme <- load10X(data_dir, features = cogaps_result, method = "CoGAPS") sme ## ----------------------------------------------------------------------------- data("curated_genes") good_gene_threshold <- 3 keep_genes <- rownames(sme)[ apply(assay(sme, "logcounts"), 1, function(x) sum(x > 0) >= good_gene_threshold) ] keep_genes <- intersect(keep_genes, curated_genes) sme <- sme[keep_genes, ] sme ## ----------------------------------------------------------------------------- spatial_params(sme) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme, feature_col = "Pattern_1", title = "Pattern_1 (immune)") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme, feature_col = "Pattern_2", title = "Pattern_2") ## ----message=FALSE, warning=FALSE--------------------------------------------- sme_ud <- find_all_hotspots(sme) ## ----------------------------------------------------------------------------- # Fraction of spots classified as hotspots per pattern hs <- hotspots(sme_ud, "undirected") pat_cols <- setdiff(colnames(hs), c("barcode", "y", "x")) vapply(pat_cols, function(p) mean(!is.na(hs[[p]])), numeric(1)) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_ud, feature_col = "Pattern_1", source = "hotspots", hotspot_type = "undirected", colors = "steelblue", title = "Pattern_1 hotspots") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_ud, feature_col = "Pattern_2", source = "hotspots", hotspot_type = "undirected", colors = "firebrick", title = "Pattern_2 hotspots") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_ud, source = "interaction", hotspot_type = "undirected", interaction_patterns = c("Pattern_1", "Pattern_2"), colors = c(Pattern_1 = "steelblue", interacting = "purple", Pattern_2 = "firebrick")) ## ----------------------------------------------------------------------------- sme_ud <- calculate_overlap_undirected(sme_ud) overlap_scores(sme_ud) analysis_type(sme_ud) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_overlap_scores(sme_ud) ## ----message=FALSE, warning=FALSE--------------------------------------------- sme_ud <- get_pairwise_interacting_genes( sme_ud, mode = "DE", analysis = "enrichment", minOverlap = 10, workers = 1, pattern_pairs = rbind(c("Pattern_1", "Pattern_2")) ) ## ----------------------------------------------------------------------------- # Detailed per-gene statistics for the Pattern_1 x Pattern_2 pair names(interactions(sme_ud)) pair1 <- interactions(sme_ud)[[1]] if (length(pair1$interacting_genes) > 0) { head(pair1$interacting_genes[[1]]) } ## ----------------------------------------------------------------------------- sme_ud <- get_im_scores(sme_ud) head(undirected_scores(sme_ud)) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ pair_col <- setdiff(colnames(undirected_scores(sme_ud)), "Gene")[1] plot_im_scores(sme_ud, interaction = pair_col, nGenes = 10) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ ims <- undirected_scores(sme_ud) top_gene <- ims[order(ims[[pair_col]], decreasing = TRUE), "Gene"][1] plot_spatial(sme_ud, feature_col = top_gene, title = sprintf("%s (top %s gene)", top_gene, pair_col)) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme, feature_col = "Pattern_1", title = "Pattern_1 (immune)") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme, feature_col = "Pattern_5", title = "Pattern_5 (cancer)") ## ----message=FALSE, warning=FALSE--------------------------------------------- sme_dir <- calculate_influence(sme) head(influence_map(sme_dir)) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, feature_col = "Pattern_1", source = "influence_map", title = "Pattern_1 influence") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, feature_col = "Pattern_5", source = "influence_map", title = "Pattern_5 influence") ## ----message=FALSE, warning=FALSE--------------------------------------------- sme_dir <- find_hotspots_gmm(sme_dir, type = "pattern") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, feature_col = "Pattern_1", source = "hotspots", hotspot_type = "pattern", colors = "steelblue", title = "Pattern_1 GMM hotspots") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, feature_col = "Pattern_5", source = "hotspots", hotspot_type = "pattern", colors = "firebrick", title = "Pattern_5 GMM hotspots") ## ----message=FALSE, warning=FALSE--------------------------------------------- sme_dir <- find_hotspots_gmm(sme_dir, type = "influence") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, feature_col = "Pattern_1", source = "hotspots", hotspot_type = "influence", colors = "steelblue", title = "Pattern_1 influence hotspots") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, feature_col = "Pattern_5", source = "hotspots", hotspot_type = "influence", colors = "firebrick", title = "Pattern_5 influence hotspots") ## ----------------------------------------------------------------------------- labs_fwd <- overlap_map(sme_dir, c("Pattern_1", "Pattern_5"), direction = "forward") table(labs_fwd, useNA = "ifany") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, source = "interaction", interaction_patterns = c("Pattern_1", "Pattern_5"), direction = "forward", colors = c(Pattern_1 = "steelblue", "Pattern_1 near Pattern_5" = "purple", "Pattern_5 influence" = "gray60")) ## ----------------------------------------------------------------------------- labs_rev <- overlap_map(sme_dir, c("Pattern_1", "Pattern_5"), direction = "reverse") table(labs_rev, useNA = "ifany") ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_spatial(sme_dir, source = "interaction", interaction_patterns = c("Pattern_1", "Pattern_5"), direction = "reverse", colors = c(Pattern_5 = "firebrick", "Pattern_5 near Pattern_1" = "goldenrod", "Pattern_1 influence" = "gray60")) ## ----------------------------------------------------------------------------- sme_dir <- calculate_overlap_directed(sme_dir) overlap_scores(sme_dir) analysis_type(sme_dir) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ plot_overlap_scores(sme_dir, title = "Directed overlap (relative abundance)") ## ----message=FALSE, warning=FALSE--------------------------------------------- sme_dir <- calculate_gene_scores_directed( sme_dir, pattern_pairs = rbind(c("Pattern_1", "Pattern_5")) ) ## ----------------------------------------------------------------------------- head(directed_scores(sme_dir)) ## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------ ds <- directed_scores(sme_dir) ci <- sort(unique(ds$cell_interaction))[1] top_gene <- ds[ds$cell_interaction == ci, ] top_gene <- top_gene[order(-top_gene$effect_size), "gene"][1] plot_spatial(sme_dir, feature_col = top_gene, title = sprintf("%s (top directed gene, %s)", top_gene, ci)) ## ----eval=FALSE--------------------------------------------------------------- # sme_ud <- sme |> # find_all_hotspots() |> # calculate_overlap_undirected() |> # get_pairwise_interacting_genes( # mode = "DE", analysis = "enrichment", # minOverlap = 10, workers = 1, # pattern_pairs = rbind(c("Pattern_1", "Pattern_2")) # ) |> # get_im_scores() # # sme_dir <- sme |> # calculate_influence() |> # find_hotspots_gmm(type = "pattern") |> # find_hotspots_gmm(type = "influence") |> # calculate_overlap_directed() |> # calculate_gene_scores_directed( # pattern_pairs = rbind(c("Pattern_1", "Pattern_5")) # ) ## ----------------------------------------------------------------------------- unlink(file.path(data_dir), recursive = TRUE) ## ----------------------------------------------------------------------------- sessionInfo()