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

## ----load_packages, message = FALSE, warning = FALSE--------------------------
library(MultiAssaySpatialExperiment)
library(SummarizedExperiment)
library(SingleCellExperiment)
library(S4Vectors)
library(sf)

## ----minimal_mase-------------------------------------------------------------
# Create a simple expression matrix
counts <- matrix(rpois(50, lambda = 10), nrow = 10, ncol = 5,
                dimnames = list(paste0("Gene", 1:10),
                               paste0("Cell", 1:5)))

# Specimen metadata
specimens <- DataFrame(
  patient_id = c("P1", "P1", "P1", "P2", "P2"),
  tissue = c("cortex", "cortex", "cortex", "medulla", "medulla"),
  row.names = paste0("Cell", 1:5)
)

# Sample map: link assay columns to specimens
sample_map <- DataFrame(
  assay = factor("rna", "rna"),
  primary = paste0("Cell", 1:5),
  colname = paste0("Cell", 1:5)
)

# Construct MASE
mase_minimal <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = counts),
  colData = specimens,
  sampleMap = sample_map
)

mase_minimal

## ----add_points_with_map------------------------------------------------------
# Suppose we have two assays: RNA and protein
rna_counts <- matrix(rpois(50, 10), nrow = 10, ncol = 5,
                    dimnames = list(paste0("Gene", 1:10), paste0("Cell", 1:5)))
prot_counts <- matrix(rpois(15, 5), nrow = 3, ncol = 5,
                     dimnames = list(paste0("Protein", 1:3), paste0("Cell", 1:5)))

# Both assays share the same specimens
specimens2 <- DataFrame(
  patient_id = c("P1", "P1", "P1", "P2", "P2"),
  row.names = paste0("Cell", 1:5)
)

# Sample map for two assays
sample_map2 <- DataFrame(
  assay = factor(rep(c("rna", "protein"), each = 5), c("rna", "protein")),
  primary = rep(paste0("Cell", 1:5), 2),
  colname = rep(paste0("Cell", 1:5), 2)
)

# Point coordinates (shared by both assays)
coords2 <- DataFrame(
  x = runif(5, 0, 100),
  y = runif(5, 0, 100),
  instance_id = paste0("Cell", 1:5)
)

# Spatial map: link assay columns to spatial points
spatial_map2 <- DataFrame(
  assay = factor(rep(c("rna", "protein"), each = 5), c("rna", "protein")),
  colname = rep(paste0("Cell", 1:5), 2),
  element_type = "points",
  region = factor(rep("coords", 10), "coords"),
  instance_id = rep(paste0("Cell", 1:5), 2)
)

# Construct MASE with spatialMap
mase_with_map <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = rna_counts, protein = prot_counts),
  colData = specimens2,
  sampleMap = sample_map2,
  points = PointsLayerList(coords = coords2),
  spatialMap = spatial_map2
)

mase_with_map

## ----add_shapes---------------------------------------------------------------
# Create simple circular boundaries around each point
centroids <- data.frame(
  x = runif(5, 0, 100),
  y = runif(5, 0, 100),
  instance_id = paste0("Cell", 1:5)
)

# Convert to sf points, then buffer to circles
centroids_sf <- st_as_sf(centroids, coords = c("x", "y"))
circles <- st_buffer(centroids_sf, dist = 5)

# Convert back to DataFrame with geometry column
boundaries <- DataFrame(
  geometry = st_geometry(circles),
  instance_id = paste0("Cell", 1:5)
)

# Add to MASE
mase_with_shapes <- mase_with_map
spatialShapes(mase_with_shapes) <- ShapesLayerList(boundaries = boundaries)

mase_with_shapes

## ----from_sce-----------------------------------------------------------------
# Create two SingleCellExperiment objects (simulating different technologies)
sce1 <- SingleCellExperiment(
  assays = list(counts = matrix(rpois(40, 10), nrow = 10, ncol = 4,
                               dimnames = list(paste0("Gene", 1:10),
                                              paste0("Cell", 1:4)))),
  colData = DataFrame(tech = "rna", cell_id = paste0("Cell", 1:4))
)

sce2 <- SingleCellExperiment(
  assays = list(counts = matrix(rpois(12, 5), nrow = 3, ncol = 4,
                               dimnames = list(paste0("Protein", 1:3),
                                              paste0("Cell", 1:4)))),
  colData = DataFrame(tech = "protein", cell_id = paste0("Cell", 1:4))
)

# Specimen metadata (shared)
specimens_sce <- DataFrame(
  patient = rep("P1", 4),
  tissue = rep("cortex", 4),
  row.names = paste0("Cell", 1:4)
)

# Sample maps
samplemap_sce <- DataFrame(
  assay = factor(rep(c("rna", "protein"), each = 4), c("rna", "protein")),
  primary = rep(paste0("Cell", 1:4), 2),
  colname = rep(paste0("Cell", 1:4), 2)
)

# Construct MASE
mase_from_sce <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = sce1, protein = sce2),
  colData = specimens_sce,
  sampleMap = samplemap_sce
)

mase_from_sce

## ----prep_mase, eval = FALSE--------------------------------------------------
# spmap <- buildSpatialMap(sample_map, region = "cells", element_type = "shapes")
# prepared <- prepMASE(
#   ExperimentList(rna = counts, protein = protein_counts),
#   colData = specimens,
#   sample_map,
#   shapes = ShapesLayerList(cells = cell_boundaries),
#   spatialMap = spmap
# )
# mase <- do.call(MultiAssaySpatialExperiment, prepared)
# drops(prepared$metadata$drops)  # rows removed during harmonization, if any

## ----check_validity-----------------------------------------------------------
validObject(mase_with_map)

## ----validation_harmonize-----------------------------------------------------
bad_sample_map <- sample_map2
bad_sample_map$colname[1] <- "NotARealColumn"
mase_dropped <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = rna_counts, protein = prot_counts),
  colData = specimens2,
  sampleMap = bad_sample_map
)

## ----validation_error, error = TRUE-------------------------------------------
try({
bad_spatial_map <- spatial_map2
bad_spatial_map$instance_id[1] <- "NotACell"
MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = rna_counts, protein = prot_counts),
  colData = specimens2,
  sampleMap = sample_map2,
  points = PointsLayerList(coords = coords2),
  spatialMap = bad_spatial_map
)
})

## ----basic-subset-setup-------------------------------------------------------
make_assay <- function() {
  SummarizedExperiment(
    assays = list(counts = matrix(rnorm(20), nrow = 5, ncol = 4,
      dimnames = list(paste0("G", 1:5), paste0("S", 1:4)))),
    # "zone", not "region": this is an arbitrary per-column covariate on the
    # assay's own colData, unrelated to spatialMap's "region" column below,
    # which names a points/shapes layer
    colData = DataFrame(zone = rep(c("core", "margin"), length.out = 4)))
}

pts <- DataFrame(x = 1:4, y = 1:4, instance_id = paste0("S", 1:4))

mase <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = make_assay(), protein = make_assay()),
  colData = DataFrame(row.names = c("P1", "P2")),
  sampleMap = DataFrame(
    assay = factor(rep(c("rna", "protein"), each = 4)),
    primary = rep(c("P1", "P1", "P2", "P2"), 2),
    colname = rep(paste0("S", 1:4), 2)
  ),
  points = PointsLayerList(coords = pts),
  spatialMap = DataFrame(
    assay = factor(rep(c("rna", "protein"), each = 4)),
    colname = rep(paste0("S", 1:4), 2),
    element_type = "points",
    region = "coords",
    instance_id = rep(paste0("S", 1:4), 2)
  )
)

vapply(experiments(mase), ncol, integer(1))

## ----subsetByColData----------------------------------------------------------
by_specimen <- mase[, "P1"]
vapply(experiments(by_specimen), ncol, integer(1))

## ----subsetByColumn-----------------------------------------------------------
by_column <- mase[, list(rna = c("S1", "S3"))]
vapply(experiments(by_column), ncol, integer(1))

## ----subsetByColumnColData----------------------------------------------------
cdf <- colData(experiments(mase)[["rna"]])
mase_core <- subsetByColumn(mase, list(rna = cdf$zone == "core"))
ncol(experiments(mase_core)[["rna"]])
nrow(spatialPoints(mase_core)[["coords"]])

## ----bracket------------------------------------------------------------------
mase[1:3, "P1", "rna", drop = FALSE]

## ----subsetByBoundingBox------------------------------------------------------
m2 <- subsetByBoundingBox(mase, xmin = 1.5, xmax = 4.5, ymin = 1.5, ymax = 4.5)
m2
spatialPoints(m2)[["coords"]]
ncol(experiments(m2)[["assay1"]])

## ----subsetByPolygon----------------------------------------------------------
poly <- st_polygon(list(matrix(
  c(1.5, 1.5, 4.5, 1.5, 4.5, 4.5, 1.5, 4.5, 1.5, 1.5),
  ncol = 2, byrow = TRUE)))
m3 <- subsetByPolygon(mase, poly)
nrow(spatialPoints(m3)[["coords"]])
ncol(experiments(m3)[["assay1"]])

## ----multiregion--------------------------------------------------------------
# Define two regions of interest
region1 <- st_polygon(list(matrix(c(1, 1, 2, 1, 2, 2, 1, 2, 1, 1), ncol = 2, byrow = TRUE)))
region2 <- st_polygon(list(matrix(c(4, 4, 5, 4, 5, 5, 4, 5, 4, 4), ncol = 2, byrow = TRUE)))

# Subset by each region
m_r1 <- subsetByPolygon(mase, region1)
m_r2 <- subsetByPolygon(mase, region2)

c(region1 = ncol(experiments(m_r1)[[1]]),
  region2 = ncol(experiments(m_r2)[[1]]))

## ----buffer_subset------------------------------------------------------------
# Reference point
ref_point <- st_point(c(3, 3))

# Create buffer (e.g., 1.5 units radius)
buffer_region <- st_buffer(ref_point, dist = 1.5)

# Subset by buffer
m_buffer <- subsetByPolygon(mase, buffer_region)
m_buffer

## ----setup-annotation---------------------------------------------------------
# Four spots at coordinates
pts_agg <- DataFrame(
  x = c(1.5, 2.5, 2.5, 3.5),
  y = c(1.5, 1.5, 2.5, 2.5),
  instance_id = paste0("S", 1:4))

# Three cell polygons
cell1 <- st_polygon(list(matrix(c(1, 1, 2, 1, 2, 2, 1, 2, 1, 1), ncol = 2, byrow = TRUE)))
cell2 <- st_polygon(list(matrix(c(2, 1, 3, 1, 3, 2, 2, 2, 2, 1), ncol = 2, byrow = TRUE)))
cell3 <- st_polygon(list(matrix(c(2, 2, 3, 2, 3, 3, 2, 3, 2, 2), ncol = 2, byrow = TRUE)))

shp_df <- DataFrame(
  instance_id = c("cell1", "cell2", "cell3"),
  geometry = st_sfc(cell1, cell2, cell3))

# Gene expression assay (3 genes × 4 spots)
expr <- matrix(c(10, 20, 5,  15,
                 30, 10, 25, 5,
                 5,  15, 20, 30),
  nrow = 3, ncol = 4,
  dimnames = list(paste0("Gene", 1:3), paste0("S", 1:4)))

# Construct MASE
mase_agg <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = expr),
  colData = DataFrame(row.names = "P1"),
  sampleMap = DataFrame(
    assay = factor("rna"),
    primary = rep("P1", 4),
    colname = paste0("S", 1:4)
  ),
  points = PointsLayerList(centroids = pts_agg),
  shapes = ShapesLayerList(cells = shp_df),
  spatialMap = DataFrame(
    assay = factor("rna"),
    colname = paste0("S", 1:4),
    element_type = "points",
    region = "centroids",
    instance_id = paste0("S", 1:4)
  )
)
mase_agg

## ----annotate-----------------------------------------------------------------
mase_agg_orig <- mase_agg
mase_agg <- annotateWithRegions(mase_agg, points = "centroids", shapes = "cells")
spatialMap(mase_agg)

## ----aggregate-count----------------------------------------------------------
agg_count <- aggregateByRegion(mase_agg, by = "cells", FUN = "count")
agg_count

## ----aggregate-sum------------------------------------------------------------
agg_sum <- aggregateByRegion(mase_agg, by = "cells", FUN = "sum")
agg_sum[["rna"]]

## ----aggregate-mean-----------------------------------------------------------
agg_mean <- aggregateByRegion(mase_agg, by = "cells", FUN = "mean")
agg_mean[["rna"]]

## ----compare_strategies-------------------------------------------------------
# Compare sum vs mean for Gene1
cells <- colnames(agg_sum[["rna"]])
data.frame(
  cell = cells,
  sum = agg_sum[["rna"]]["Gene1", ],
  mean = agg_mean[["rna"]]["Gene1", ],
  n_spots = agg_count$count[match(cells, agg_count$cells)]
)

## ----vis_aggregated, fig.width = 8, fig.height = 4----------------------------
# Get aggregated data for Gene1
gene1_expr <- agg_sum[["rna"]]["Gene1", ]

# Get cell polygons as sf
cells_sf <- st_sf(
  instance_id = shp_df$instance_id,
  geometry = shp_df$geometry
)

# Add expression values
cells_sf$Gene1_sum <- gene1_expr[cells_sf$instance_id]

# Plot
par(mfrow = c(1, 2))

# Original points
plot(pts_agg$x, pts_agg$y, 
     pch = 16, cex = 2,
     col = rainbow(4)[rank(expr["Gene1", ])],
     main = "Gene1: Original spots",
     xlab = "x", ylab = "y")
text(pts_agg$x, pts_agg$y, paste0("S", 1:4), pos = 3, cex = 0.7)

# Aggregated cells
plot(st_geometry(cells_sf), 
     col = rainbow(3)[rank(cells_sf$Gene1_sum)],
     main = "Gene1: Aggregated by cell",
     xlab = "x", ylab = "y")
text(st_coordinates(st_centroid(cells_sf)), 
     labels = round(cells_sf$Gene1_sum, 1),
     cex = 0.8)

## ----rasterize----------------------------------------------------------------
# Example: 3 cell polygons
cell1_rast <- st_polygon(list(matrix(c(0, 0, 1, 0, 1, 1, 0, 1, 0, 0), ncol = 2, byrow = TRUE)))
cell2_rast <- st_polygon(list(matrix(c(1, 0, 2, 0, 2, 1, 1, 1, 1, 0), ncol = 2, byrow = TRUE)))
cell3_rast <- st_polygon(list(matrix(c(0, 1, 1, 1, 1, 2, 0, 2, 0, 1), ncol = 2, byrow = TRUE)))

shapes_df <- DataFrame(
  instance_id = c("cell1", "cell2", "cell3"),
  geometry = st_sfc(cell1_rast, cell2_rast, cell3_rast)
)

## ----vectorize----------------------------------------------------------------
# Example label matrix (3 cells)
label_matrix <- matrix(c(
  1, 1, 2, 2,
  1, 1, 2, 2,
  3, 3, 0, 0,
  3, 3, 0, 0
), nrow = 4, byrow = TRUE)

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