## ---- 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)

## -----------------------------------------------------------------------------
# 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))

## -----------------------------------------------------------------------------
# Example from SpatialExperiment
dir <- system.file(
  file.path("extdata", "10xVisium"), 
  package = "SpatialExperiment")
  
sample_ids <- c("section1", "section2")
samples <- file.path(dir, sample_ids, "outs")
  
list.files(samples[1])
list.files(file.path(samples[1], "spatial"))
file.path(samples[1], "raw_feature_bc_matrix")

## -----------------------------------------------------------------------------
(sfe3 <- read10xVisiumSFE(samples, sample_ids, type = "sparse", data = "raw",
                         load = FALSE))

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

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

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

## -----------------------------------------------------------------------------
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")

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