## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#  BiocManager::install("SpatialFeatureExperiment")

## ----echo=FALSE, out.width = "100%", fig.cap="Schematics of the SFE object", fig.alt="SpatialFeatureExperiment expands on SpatialExperiment by adding column, row, and annotation geometries and spatial graphs. This is explained in detail in the following paragraphs."----
knitr::include_graphics("sfe_schematics.png")

## ----setup--------------------------------------------------------------------
library(SpatialFeatureExperiment)
library(SpatialExperiment)
library(SFEData)
library(sf)
library(Matrix)
library(RBioFormats)

## -----------------------------------------------------------------------------
# Example dataset
(sfe <- McKellarMuscleData(dataset = "small"))

## -----------------------------------------------------------------------------
# Get Visium spot polygons
(spots <- dimGeometry(sfe, "spotPoly", MARGIN = 2))

## -----------------------------------------------------------------------------
plot(st_geometry(spots))

## -----------------------------------------------------------------------------
# Setter
dimGeometry(sfe, "foobar", MARGIN = 2) <- spots

## -----------------------------------------------------------------------------
# Getter, all geometries of one margin
(cgs <- dimGeometries(sfe, MARGIN = 2))

## -----------------------------------------------------------------------------
# Setter, all geometries
dimGeometries(sfe, MARGIN = 2) <- cgs

## -----------------------------------------------------------------------------
(cg_names <- dimGeometryNames(sfe, MARGIN = 2))

## -----------------------------------------------------------------------------
# Setter
dimGeometryNames(sfe, MARGIN = 2) <- cg_names

## -----------------------------------------------------------------------------
# Getter
(spots <- spotPoly(sfe))

## -----------------------------------------------------------------------------
# Setter
spotPoly(sfe) <- spots

## -----------------------------------------------------------------------------
# Getter, by name or index
(tb <- annotGeometry(sfe, "tissueBoundary"))

## -----------------------------------------------------------------------------
plot(st_geometry(tb))

## -----------------------------------------------------------------------------
# Setter, by name or index
annotGeometry(sfe, "tissueBoundary") <- tb

## -----------------------------------------------------------------------------
# Get all annoation geometries as named list
ags <- annotGeometries(sfe)

## -----------------------------------------------------------------------------
# Set all annotation geometries with a named list
annotGeometries(sfe) <- ags

## -----------------------------------------------------------------------------
# Get names of annotation geometries
(ag_names <- annotGeometryNames(sfe))

## -----------------------------------------------------------------------------
# Set names
annotGeometryNames(sfe) <- ag_names

## -----------------------------------------------------------------------------
# Getter
(tb <- tissueBoundary(sfe))

## -----------------------------------------------------------------------------
# Setter
tissueBoundary(sfe) <- tb

## -----------------------------------------------------------------------------
(g <- findSpatialNeighbors(sfe, MARGIN = 2, method = "tri2nb"))

## -----------------------------------------------------------------------------
plot(g, coords = spatialCoords(sfe))

## -----------------------------------------------------------------------------
# Set graph by name
spatialGraph(sfe, "graph1", MARGIN = 2) <- g
# Or equivalently
colGraph(sfe, "graph1") <- g

## -----------------------------------------------------------------------------
# Get graph by name
g <- spatialGraph(sfe, "graph1", MARGIN = 2L)
# Or equivalently
g <- colGraph(sfe, "graph1")
g

## -----------------------------------------------------------------------------
colGraph(sfe, "visium") <- findVisiumGraph(sfe)

## -----------------------------------------------------------------------------
plot(colGraph(sfe, "visium"), coords = spatialCoords(sfe))

## -----------------------------------------------------------------------------
colGraphs(sfe)

## -----------------------------------------------------------------------------
colGraphNames(sfe)

## -----------------------------------------------------------------------------
# Construct toy dataset with 2 samples
sfe1 <- McKellarMuscleData(dataset = "small")
sfe2 <- McKellarMuscleData(dataset = "small2")
spotPoly(sfe2)$sample_id <- "sample02"
(sfe_combined <- cbind(sfe1, sfe2))

## -----------------------------------------------------------------------------
sampleIDs(sfe_combined)

## -----------------------------------------------------------------------------
# Only get the geometries for the second sample
(spots2 <- colGeometry(sfe_combined, "spotPoly", sample_id = "sample02"))

## -----------------------------------------------------------------------------
# Only set the geometries for the second sample
# Leaving geometries of the first sample intact
colGeometry(sfe_combined, "spotPoly", sample_id = "sample02") <- spots2

## -----------------------------------------------------------------------------
# Set graph only for the second sample
colGraph(sfe_combined, "foo", sample_id = "sample02") <- 
  findSpatialNeighbors(sfe_combined, sample_id = "sample02")

## -----------------------------------------------------------------------------
# Get graph only for the second sample
colGraph(sfe_combined, "foo", sample_id = "sample02")

## -----------------------------------------------------------------------------
# Set graph of the same name for both samples
# The graphs are computed separately for each sample
colGraphs(sfe_combined, sample_id = "all", name = "visium") <- 
  findVisiumGraph(sfe_combined, sample_id = "all")

## -----------------------------------------------------------------------------
# Get multiple graphs of the same name
colGraphs(sfe_combined, sample_id = "all", name = "visium")

## -----------------------------------------------------------------------------
# Or just all graphs of the margin
colGraphs(sfe_combined, sample_id = "all")

## -----------------------------------------------------------------------------
sfe_combined <- changeSampleIDs(sfe, replacement = c(Vis5A = "foo", sample02 = "bar"))
sfe_combined

## -----------------------------------------------------------------------------
# Visium barcode location from Space Ranger
data("visium_row_col")
coords1 <- visium_row_col[visium_row_col$col < 6 & visium_row_col$row < 6,]
coords1$row <- coords1$row * sqrt(3)

# Random toy sparse matrix
set.seed(29)
col_inds <- sample(1:13, 13)
row_inds <- sample(1:5, 13, replace = TRUE)
values <- sample(1:5, 13, replace = TRUE)
mat <- sparseMatrix(i = row_inds, j = col_inds, x = values)
colnames(mat) <- coords1$barcode
rownames(mat) <- sample(LETTERS, 5)

## -----------------------------------------------------------------------------
sfe3 <- SpatialFeatureExperiment(list(counts = mat), colData = coords1,
                                spatialCoordsNames = c("col", "row"),
                                spotDiameter = 0.7)

## -----------------------------------------------------------------------------
# Convert regular data frame with coordinates to sf data frame
cg <- df2sf(coords1[,c("col", "row")], c("col", "row"), spotDiameter = 0.7)
rownames(cg) <- colnames(mat)
sfe3 <- SpatialFeatureExperiment(list(counts = mat), colGeometries = list(foo = cg))

## -----------------------------------------------------------------------------
dir <- system.file("extdata", package = "SpatialFeatureExperiment")
sample_ids <- c("sample01", "sample02")
samples <- file.path(dir, sample_ids)

## -----------------------------------------------------------------------------
list.files(file.path(samples[1], "outs"))

## -----------------------------------------------------------------------------
list.files(file.path(samples[1], "outs", "filtered_feature_bc_matrix"))

## -----------------------------------------------------------------------------
list.files(file.path(samples[1], "outs", "spatial"))

## -----------------------------------------------------------------------------
(sfe3 <- read10xVisiumSFE(samples, sample_id = sample_ids, type = "sparse", 
                          data = "filtered", images = "hires"))

## -----------------------------------------------------------------------------
(sfe3 <- read10xVisiumSFE(samples, sample_id = sample_ids, type = "sparse", 
                          data = "filtered", images = "hires", unit = "micron"))

## -----------------------------------------------------------------------------
unit(sfe3)

## -----------------------------------------------------------------------------
class(imgRaster(getImg(sfe3)))

## -----------------------------------------------------------------------------
fp <- tempdir()
dir_use <- VizgenOutput(file_path = file.path(fp, "vizgen"))
list.files(dir_use)

## -----------------------------------------------------------------------------
(sfe_mer <- readVizgen(dir_use, z = 3L, image = "PolyT", add_molecules = TRUE))

## -----------------------------------------------------------------------------
list.files(dir_use)

## -----------------------------------------------------------------------------
dir_use <- XeniumOutput("v2", file_path = file.path(fp, "xenium"))
list.files(dir_use)

## -----------------------------------------------------------------------------
# RBioFormats issue: https://github.com/aoles/RBioFormats/issues/42
try(sfe_xen <- readXenium(dir_use, add_molecules = TRUE))
(sfe_xen <- readXenium(dir_use, add_molecules = TRUE))

## -----------------------------------------------------------------------------
list.files(dir_use)

## -----------------------------------------------------------------------------
dir_use <- CosMXOutput(file_path = file.path(fp, "cosmx"))
list.files(dir_use)

## -----------------------------------------------------------------------------
(sfe_cosmx <- readCosMX(dir_use, add_molecules = TRUE))

## -----------------------------------------------------------------------------
list.files(dir_use)

## -----------------------------------------------------------------------------
spe <- read10xVisium(samples, sample_ids, type = "sparse", data = "filtered", 
  images = "hires", load = FALSE)

## -----------------------------------------------------------------------------
colnames(spe) <- make.unique(colnames(spe), sep = "-")
rownames(spatialCoords(spe)) <- colnames(spe)

## -----------------------------------------------------------------------------
(sfe3 <- toSpatialFeatureExperiment(spe))

## -----------------------------------------------------------------------------
dir_extdata <- system.file("extdata", package = "SpatialFeatureExperiment")
obj_vis <- readRDS(file.path(dir_extdata, "seu_vis_toy.rds"))

## -----------------------------------------------------------------------------
sfe_conv_vis <-
  toSpatialFeatureExperiment(x = obj_vis, 
                             image_scalefactors = "lowres",
                             unit = "micron",
                             BPPARAM = BPPARAM)
sfe_conv_vis

## -----------------------------------------------------------------------------
sfe_combined <- cbind(sfe1, sfe2)

## -----------------------------------------------------------------------------
(sfe_subset <- sfe[1:10, 1:10, drop = TRUE])

## ----eval=FALSE---------------------------------------------------------------
#  # Will give warning because graph reconstruction fails
#  sfe_subset <- sfe[1:10, 1:10]

## -----------------------------------------------------------------------------
# Before
plot(st_geometry(tissueBoundary(sfe)))
plot(spotPoly(sfe), col = "gray", add = TRUE)

## -----------------------------------------------------------------------------
sfe_in_tissue <- crop(sfe, y = tissueBoundary(sfe), colGeometryName = "spotPoly")

## -----------------------------------------------------------------------------
# After
plot(st_geometry(tissueBoundary(sfe)))
plot(spotPoly(sfe_in_tissue), col = "gray", add = TRUE)

## -----------------------------------------------------------------------------
sfe_cropped <- crop(sfe, colGeometryName = "spotPoly", sample_id = "Vis5A",
                    xmin = 5500, xmax = 6500, ymin = 13500, ymax = 14500)

## -----------------------------------------------------------------------------
# Get logical vector
colData(sfe)$in_tissue <- annotPred(sfe, colGeometryName = "spotPoly", 
                                    annotGeometryName = "tissueBoundary",
                                    sample_id = "Vis5A")
# Get the number of nuclei per Visium spot
colData(sfe)$n_nuclei <- annotNPred(sfe, "spotPoly", annotGeometryName = "nuclei")
# Get geometries of intersections of Visium spots and myofibers
spot_intersections <- annotOp(sfe, colGeometryName = "spotPoly", 
                              annotGeometryName = "myofiber_simplified")

## -----------------------------------------------------------------------------
bbox(sfe, sample_id = "Vis5A")

## -----------------------------------------------------------------------------
sfe_moved <- removeEmptySpace(sfe, sample_id = "Vis5A")

## ----fig.width=1, fig.height=1------------------------------------------------
img <- getImg(sfe3, image_id = "hires")
plot(imgRaster(img))
plot(transposeImg(img) |> imgRaster())

## ----fig.width=1, fig.height=1------------------------------------------------
plot(mirrorImg(img, direction = "vertical") |> imgRaster())
plot(mirrorImg(img, direction = "horizontal") |> imgRaster())

## -----------------------------------------------------------------------------
sfe3 <- mirrorImg(sfe3, sample_id = "sample01", image_id = "hires")

## -----------------------------------------------------------------------------
sfe_mirrored <- mirror(sfe_in_tissue)
sfe_transposed <- transpose(sfe_in_tissue)

## ----fig.width=6, fig.height=2------------------------------------------------
par(mfrow = c(1, 3), mar = rep(1.5, 4))
plot(st_geometry(tissueBoundary(sfe_in_tissue)))
plot(spotPoly(sfe_in_tissue), col = "gray", add = TRUE)

plot(st_geometry(tissueBoundary(sfe_mirrored)))
plot(spotPoly(sfe_mirrored), col = "gray", add = TRUE)

plot(st_geometry(tissueBoundary(sfe_transposed)))
plot(spotPoly(sfe_transposed), col = "gray", add = TRUE)

## -----------------------------------------------------------------------------
# Clean up
unlink(file.path(fp, "vizgen"), recursive = TRUE)
unlink(file.path(fp, "xenium"), recursive = TRUE)
unlink(file.path(fp, "cosmx"), recursive = TRUE)

## -----------------------------------------------------------------------------
sessionInfo()