## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  cache = TRUE,
  fig.width = 7,
  fig.height = 5
)

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

# For visualization
library(ggplot2)
set.seed(123)

## ----read_xenium, eval = FALSE------------------------------------------------
# mase_xenium <- readXeniumMASE(
#   data_dir = "path/to/xenium_output",
#   sample_id = "sample1",
#   segmentations = "cell",
#   add_transcripts = FALSE  # Skip transcripts for faster loading
# )

## ----read_visium, eval = FALSE------------------------------------------------
# mase_visium <- readVisiumMASE(
#   data_dir = "path/to/spaceranger_output",
#   sample_id = "sample1",
#   images = TRUE,
#   data = "filtered"
# )

## ----read_visium_hd, eval = FALSE---------------------------------------------
# # Read Visium HD at 8 µm resolution
# mase_visium_hd <- readVisiumHDMASE(
#   data_dir = "path/to/spaceranger_hd_output",
#   sample_id = "sample_hd",
#   bin_size = "008",
#   images = TRUE
# )

## ----read_cosmx, eval = FALSE-------------------------------------------------
# mase_cosmx <- readCosMxMASE(
#   data_dir = "path/to/cosmx_output",
#   sample_id = "cosmx_sample",
#   fov_ids = NULL  # NULL loads all FOVs; or c("1", "3", "5")
# )

## ----read_merfish, eval = FALSE-----------------------------------------------
# mase_merscope <- readMERSCOPEMASE(
#   data_dir = "path/to/vizgen_output",
#   sample_id = "merscope_sample",
#   segmentation = "cellpose",
#   load_transcripts = FALSE
# )

## ----from_spe, eval = FALSE---------------------------------------------------
# library(SpatialExperiment)
# 
# # spe <- read10xVisium("path/to/visium")
# mase_from_spe <- as(spe, "MultiAssaySpatialExperiment")

## ----from_sfe, eval = FALSE---------------------------------------------------
# library(SpatialFeatureExperiment)
# 
# # sfe <- read10xVisiumSFE("path/to/visium")
# mase_from_sfe <- as(sfe, "MultiAssaySpatialExperiment")

## ----multi_sample_read, eval = FALSE------------------------------------------
# mase1 <- readXeniumMASE(data_dir = "path/to/sample1", sample_id = "P001_xenium")
# mase2 <- readXeniumMASE(data_dir = "path/to/sample2", sample_id = "P002_xenium")
# 
# mase_multi <- c(mase1, mase2)

## ----prepare_xenium-----------------------------------------------------------
# Xenium data: single-cell resolution
# Simulate 500 cells, 50 genes (subset of full panel)
n_cells <- 500
n_genes_xenium <- 50

xenium_counts <- matrix(
  rpois(n_cells * n_genes_xenium, lambda = 5),
  nrow = n_genes_xenium,
  ncol = n_cells,
  dimnames = list(
    paste0("Gene", 1:n_genes_xenium),
    paste0("Cell", 1:n_cells)
  )
)

# Cell centroids (tissue coordinates in micrometers)
xenium_coords <- DataFrame(
  x = runif(n_cells, 0, 1000),
  y = runif(n_cells, 0, 1000),
  instance_id = paste0("Cell", 1:n_cells)
)

# Cell type annotations (from clustering or known markers)
cell_types <- sample(
  c("Neuron", "Astrocyte", "Microglia", "Oligodendrocyte"),
  n_cells,
  replace = TRUE,
  prob = c(0.5, 0.2, 0.15, 0.15)
)

# Create SingleCellExperiment
xenium_sce <- SingleCellExperiment(
  assays = list(counts = xenium_counts),
  colData = DataFrame(
    cell_type = cell_types,
    tech = "xenium"
  )
)

## ----prepare_visium-----------------------------------------------------------
# Visium data: spot-level resolution
# Simulate 100 spots, 1000 genes (whole transcriptome, downsampled)
n_spots <- 100
n_genes_visium <- 1000

# Create gene names: first 50 overlap with Xenium, rest are Visium-only
visium_gene_names <- c(
  paste0("Gene", 1:n_genes_xenium),  # Overlapping genes (50)
  paste0("GeneV", 1:(n_genes_visium - n_genes_xenium))  # Visium-only (950)
)

visium_counts <- matrix(
  rpois(n_spots * n_genes_visium, lambda = 10),
  nrow = n_genes_visium,
  ncol = n_spots,
  dimnames = list(
    visium_gene_names,
    paste0("Spot", 1:n_spots)
  )
)

# Spot centroids (same coordinate system as Xenium)
visium_coords <- DataFrame(
  x = runif(n_spots, 0, 1000),
  y = runif(n_spots, 0, 1000),
  instance_id = paste0("Spot", 1:n_spots)
)

# Create SummarizedExperiment
visium_se <- SummarizedExperiment(
  assays = list(counts = visium_counts),
  colData = DataFrame(
    tech = rep("visium", n_spots),
    row.names = paste0("Spot", 1:n_spots)
  )
)

## ----build_mase---------------------------------------------------------------
# Specimen metadata: one row per biological unit (here, one tissue section)
specimen_id <- "section_01"
specimens <- DataFrame(
  patient_id = "P001",
  tissue = "mouse_brain",
  section_id = specimen_id,
  row.names = specimen_id
)

n_total <- n_cells + n_spots

# Sample map: link each assay column (cell or spot) to the shared specimen
sample_map <- DataFrame(
  assay = factor(
    c(rep("xenium", n_cells), rep("visium", n_spots)),
    c("xenium", "visium")
  ),
  primary = rep(specimen_id, n_total),
  colname = c(paste0("Cell", 1:n_cells), paste0("Spot", 1:n_spots))
)

# Spatial map: link assay columns to spatial coordinates
spatial_map <- DataFrame(
  assay = factor(
    c(rep("xenium", n_cells), rep("visium", n_spots)),
    c("xenium", "visium")
  ),
  colname = c(paste0("Cell", 1:n_cells), paste0("Spot", 1:n_spots)),
  element_type = rep("points", n_total),
  region = factor(
    c(rep("xenium_coords", n_cells), rep("visium_coords", n_spots)),
    c("xenium_coords", "visium_coords")
  ),
  instance_id = c(paste0("Cell", 1:n_cells), paste0("Spot", 1:n_spots))
)

# Construct MASE
mase <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(
    xenium = xenium_sce,
    visium = visium_se
  ),
  colData = specimens,
  sampleMap = sample_map,
  points = PointsLayerList(
    xenium_coords = xenium_coords,
    visium_coords = visium_coords
  ),
  spatialMap = spatial_map
)

mase

## ----vis_spatial, fig.width = 10, fig.height = 5------------------------------
# Extract coordinates
xenium_pts <- spatialPoint(mase, "xenium_coords")
visium_pts <- spatialPoint(mase, "visium_coords")

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

# Xenium cells colored by type
# widen the x range so the legend sits beside the points rather than on top of them
plot(xenium_pts$x, xenium_pts$y,
     col = as.factor(cell_types),
     pch = 16, cex = 0.5,
     main = "Xenium: Single cells",
     xlab = "x (µm)", ylab = "y (µm)",
     xlim = c(min(xenium_pts$x), max(xenium_pts$x) * 1.35))
legend("topright", legend = levels(as.factor(cell_types)),
       col = seq_along(levels(as.factor(cell_types))), pch = 16,
       cex = 0.7, bg = "white", box.col = "grey70")

# Visium spots
plot(visium_pts$x, visium_pts$y,
     pch = 21, bg = "lightblue", cex = 2,
     main = "Visium: Spots (55 µm)",
     xlab = "x (µm)", ylab = "y (µm)")

## ----gene_overlap-------------------------------------------------------------
xenium_genes <- rownames(experiments(mase)$xenium)
visium_genes <- rownames(experiments(mase)$visium)

shared_genes <- intersect(xenium_genes, visium_genes)

c(shared = length(shared_genes),
  xenium_only = length(setdiff(xenium_genes, visium_genes)),
  visium_only = length(setdiff(visium_genes, xenium_genes)))

## ----spatial_annotation-------------------------------------------------------
# Create circular regions for Visium spots (55 µm diameter = 27.5 µm radius)
visium_circles <- st_buffer(
  st_sf(visium_pts,
        geometry = st_as_sfc(paste0("POINT(", visium_pts$x, " ", visium_pts$y, ")"))),
  dist = 27.5
)

visium_shapes <- DataFrame(
  geometry = st_geometry(visium_circles),
  instance_id = visium_pts$instance_id
)

# Add shapes to MASE
spatialShapes(mase) <- ShapesLayerList(visium_spots = visium_shapes)

# For each Visium spot, find which Xenium cells are inside
xenium_sf <- st_sf(xenium_pts,
                   geometry = st_as_sfc(paste0("POINT(", xenium_pts$x, " ", xenium_pts$y, ")")))
visium_sf <- st_sf(
  instance_id = visium_shapes$instance_id,
  geometry = visium_shapes$geometry
)

# Spatial intersection
intersections <- st_intersects(visium_sf, xenium_sf)

# Calculate cell type proportions per spot
spot_cell_types <- lapply(seq_along(intersections), function(i) {
  cell_indices <- intersections[[i]]
  if (length(cell_indices) == 0) {
    return(data.frame(
      spot = visium_pts$instance_id[i],
      Neuron = 0, Astrocyte = 0, Microglia = 0, Oligodendrocyte = 0,
      n_cells = 0
    ))
  }
  types <- cell_types[cell_indices]
  props <- base::table(types) / length(types)
  data.frame(
    spot = visium_pts$instance_id[i],
    Neuron = as.numeric(props["Neuron"]),
    Astrocyte = as.numeric(props["Astrocyte"]),
    Microglia = as.numeric(props["Microglia"]),
    Oligodendrocyte = as.numeric(props["Oligodendrocyte"]),
    n_cells = length(cell_indices),
    stringsAsFactors = FALSE
  )
})

spot_annotations <- do.call(rbind, spot_cell_types)
spot_annotations[is.na(spot_annotations)] <- 0

head(spot_annotations)

## ----vis_deconv, fig.width = 10, fig.height = 8-------------------------------
# Plot cell type proportions per spot
par(mfrow = c(2, 2))

for (cell_type in c("Neuron", "Astrocyte", "Microglia", "Oligodendrocyte")) {
  plot(visium_pts$x, visium_pts$y,
       pch = 21, cex = 3,
       bg = rgb(spot_annotations[[cell_type]], 0, 0, alpha = 0.7),
       main = paste(cell_type, "proportion"),
       xlab = "x (µm)", ylab = "y (µm)")
}

## ----define_regions-----------------------------------------------------------
# Define two regions: "cortex" (y > 500) and "subcortex" (y <= 500)
# In practice, these would be anatomical annotations
region_polygons <- list(
  cortex = st_polygon(list(cbind(
    c(0, 1000, 1000, 0, 0),
    c(500, 500, 1000, 1000, 500)
  ))),
  subcortex = st_polygon(list(cbind(
    c(0, 1000, 1000, 0, 0),
    c(0, 0, 500, 500, 0)
  )))
)

regions_sf <- st_sf(
  region_id = c("cortex", "subcortex"),
  geometry = st_sfc(region_polygons)
)

# Convert to DataFrame for MASE
anatomical_regions_df <- DataFrame(
  geometry = st_geometry(regions_sf),
  instance_id = c("cortex", "subcortex")
)

# Add to MASE, as a shapes layer named "anatomical_regions" (not "regions": that
# name is reserved for spatialMap's own `region` column, which names the points
# layer each row's instance_id comes from, "xenium_coords" here)
spatialShapes(mase) <- ShapesLayerList(
  visium_spots = visium_shapes,
  anatomical_regions = anatomical_regions_df
)

## ----aggregate_by_region------------------------------------------------------
# Annotate Xenium centroids with anatomical regions, then aggregate expression.
# Run on the full MASE so the anatomical_regions layer is retained
# (subsetByAssay drops shapes not referenced in spatialMap).
mase <- annotateWithRegions(mase,
    points = "xenium_coords", shapes = "anatomical_regions")
xenium_by_region <- aggregateByRegion(mase, by = "anatomical_regions",
    assays = "xenium", FUN = "sum")
xenium_by_region[["xenium"]][1:5, ]

## ----compare_assays-----------------------------------------------------------
# Extract counts for shared genes
xenium_shared <- assay(experiments(mase)$xenium, "counts")[shared_genes, ]
visium_shared <- assay(experiments(mase)$visium, "counts")[shared_genes, ]

# Calculate mean expression per gene
xenium_means <- rowMeans(xenium_shared)
visium_means <- rowMeans(visium_shared)

# Plot correlation
plot(log1p(xenium_means), log1p(visium_means),
     pch = 16, col = rgb(0, 0, 0, 0.5),
     xlab = "log1p(mean Xenium counts)",
     ylab = "log1p(mean Visium counts)",
     main = "Gene expression correlation")
abline(0, 1, col = "red", lty = 2)

# Correlation
cor_val <- cor(xenium_means, visium_means, method = "spearman")
text(0.5, max(log1p(visium_means)) * 0.9,
     paste("Spearman rho =", round(cor_val, 3)),
     pos = 4)

## ----vis_to_microns-----------------------------------------------------------
# Visium spot coordinates (grid indices from Space Ranger)
visium_grid <- DataFrame(
  row_idx = c(0, 0, 1, 1, 2, 2),
  col_idx = c(0, 1, 0, 1, 0, 1),
  instance_id = paste0("Spot", 1:6)
)

# Convert to microns (approximate)
# Spot spacing: ~100 µm horizontally, ~87 µm vertically (hexagonal)
visium_microns <- DataFrame(
  x = visium_grid$col_idx * 100,  # Column → x
  y = visium_grid$row_idx * 87,   # Row → y (hexagonal offset)
  instance_id = visium_grid$instance_id
)

visium_microns

## ----affine_transform---------------------------------------------------------
# Define affine matrix: translate by (50, 30), scale by 1.2, no rotation
# Format: 2x3 matrix [a b c; d e f] where:
#   x' = a*x + b*y + c
#   y' = d*x + e*y + f

affine_matrix <- matrix(c(
  1.2, 0.0, 50,   # x' = 1.2*x + 0*y + 50 (scale + translate)
  0.0, 1.2, 30    # y' = 0*x + 1.2*y + 30 (scale + translate)
), nrow = 2, byrow = TRUE)

# Apply transformation
apply_affine <- function(pts, mat) {
  x_new <- mat[1, 1] * pts$x + mat[1, 2] * pts$y + mat[1, 3]
  y_new <- mat[2, 1] * pts$x + mat[2, 2] * pts$y + mat[2, 3]
  DataFrame(x = x_new, y = y_new, instance_id = pts$instance_id)
}

# Simulate Xenium coordinates
xenium_original <- DataFrame(
  x = c(100, 150, 200, 250, 300),
  y = c(100, 150, 200, 150, 100),
  instance_id = paste0("Cell", 1:5)
)

# Transform Xenium to align with Visium
xenium_aligned <- apply_affine(xenium_original, affine_matrix)

xenium_aligned

## ----landmark_alignment-------------------------------------------------------
# Example: three landmarks in both Xenium and Visium
# In practice, these would be manually identified features (blood vessels, etc.)

xenium_landmarks <- matrix(c(
  100, 100,
  200, 150,
  300, 200
), ncol = 2, byrow = TRUE)

visium_landmarks <- matrix(c(
  50, 50,
  150, 100,
  250, 150
), ncol = 2, byrow = TRUE)

## ----store_transformed_coords-------------------------------------------------
# Create MASE with both coordinate systems
# Original Xenium coordinates
pts_original <- PointsLayerList(xenium_original = xenium_original)

# Aligned Xenium coordinates
pts_aligned <- PointsLayerList(xenium_aligned = xenium_aligned)

# Visium coordinates
pts_visium <- PointsLayerList(visium_coords = visium_microns)

## ----vis_alignment, fig.width = 10, fig.height = 5----------------------------
par(mfrow = c(1, 2))

# Before alignment
plot(xenium_original$x, xenium_original$y,
     pch = 16, col = "red", cex = 1.5,
     main = "Before alignment",
     xlab = "x (µm)", ylab = "y (µm)",
     xlim = c(0, 400), ylim = c(0, 350))
points(visium_microns$x, visium_microns$y,
       pch = 21, bg = "blue", cex = 3)
legend("topright", legend = c("Xenium (cells)", "Visium (spots)"),
       col = c("red", "blue"), pch = c(16, 21), pt.cex = c(1.5, 2))

# After alignment
plot(xenium_aligned$x, xenium_aligned$y,
     pch = 16, col = "red", cex = 1.5,
     main = "After alignment",
     xlab = "x (µm)", ylab = "y (µm)",
     xlim = c(0, 400), ylim = c(0, 350))
points(visium_microns$x, visium_microns$y,
       pch = 21, bg = "blue", cex = 3)
legend("topright", legend = c("Xenium (aligned)", "Visium (spots)"),
       col = c("red", "blue"), pch = c(16, 21), pt.cex = c(1.5, 2))

## ----multi_sample, eval = FALSE-----------------------------------------------
# specimens_multi <- DataFrame(
#   patient_id = c(rep("P001", 150), rep("P002", 120), rep("P003", 180)),
#   tissue = "mouse_brain",
#   condition = c(rep("control", 150), rep("treated", 120), rep("treated", 180)),
#   row.names = paste0("Sample", 1:450)
# )
# 
# # sampleMap is built exactly as in the single-sample example above, with one
# # row per assay column across all three patients.

## ----multi_sample_subset, eval = FALSE----------------------------------------
# mase[, colData(mase)$patient_id == "P001", ]   # one patient, every assay
# mase[, colData(mase)$condition == "control", ] # one arm of the study
# subsetByAssay(mase, "xenium")                  # one modality, every patient

## ----workflow_qc, eval = FALSE------------------------------------------------
# mase <- readXeniumMASE(data_dir = "path/to/xenium_output", sample_id = "sample1")
# mase_qc <- mase[, colData(mase)$qc_pass, ]

## ----workflow_roi, eval = FALSE-----------------------------------------------
# library(sf)
# 
# roi <- st_polygon(list(matrix(
#   c(100, 200, 500, 200, 500, 600, 100, 600, 100, 200),
#   ncol = 2, byrow = TRUE)))
# mase_roi <- subsetByPolygon(mase_qc, roi)

## ----workflow_agg, eval = FALSE-----------------------------------------------
# mase_annotated <- annotateWithRegions(mase_roi,
#     points = "centroids", shapes = "cells")
# 
# cell_expr <- aggregateByRegion(mase_annotated, by = "cells", FUN = "sum")

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