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

## ----installation, eval=FALSE-------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }
# BiocManager::install("panoramic")

## ----libraries, message=FALSE, warning=FALSE----------------------------------
library(panoramic)
library(SummarizedExperiment)
library(BiocParallel)
library(dplyr)
library(ggplot2)

## ----simulate_data, message=FALSE, warning=FALSE------------------------------
set.seed(123)
toy <- panoramic_simulate_dataset(
  n_group1 = 3,
  n_group2 = 3,
  n_cells_group1 = 200,
  n_cells_group2 = 350,
  group_labels = c("group1", "group2"),
  scenario_group1 = "random",
  scenario_group2 = "colocalized",
  seed = 123
)

spe_list <- toy$spe_list
design <- toy$design
design

## ----plot_samples, message=FALSE, warning=FALSE-------------------------------
plot_df <- do.call(rbind, lapply(spe_list, function(spe) {
  coords <- SpatialExperiment::spatialCoords(spe)
  data.frame(
    x = coords[, 1],
    y = coords[, 2],
    cell_type = SummarizedExperiment::colData(spe)$cell_type,
    sample_id = SummarizedExperiment::colData(spe)$sample_id,
    stringsAsFactors = FALSE
  )
}))

ggplot(plot_df, aes(x = x, y = y, color = cell_type)) +
  geom_point(size = 0.7, alpha = 0.7) +
  facet_wrap(~ sample_id, nrow = 2) +
  coord_equal() +
  theme_bw() +
  labs(x = "x", y = "y", color = "Cell type")

## ----prepare, message=FALSE, warning=FALSE------------------------------------
prep <- panoramic_prepare(
  spe_list = spe_list,
  design = design,
  cell_type = "cell_type",
  min_cells = 5,
  window = "concave"
)
names(S4Vectors::metadata(prep[[1]])$panoramic)

## ----spatial_stats, message=FALSE, warning=FALSE------------------------------
radii_um <- c(25)

se <- panoramic_spatialstats(
  prep = prep,
  radii_um = radii_um,
  stat = "local_comp_enrichment",
  nsim = 30,
  seed = 123,
  BPPARAM = BiocParallel::SerialParam()
)
dim(se)
head(as.data.frame(rowData(se))[, c("ct1", "ct2", "radius_um", "stat")])

## ----meta_mv, message=FALSE, warning=FALSE------------------------------------
# Here patient_col and sample_col are both "sample" in toy data.
# In real data, patient_col can differ from sample_col (e.g., multiple cores per patient).

se_meta <- panoramic_meta_mv(
  se,
  patient_col = "sample",
  group_col = "group",
  sample_col = "sample",
  tau_structure = "patient",
  method = "REML",
  group_tau2 = "separate"
)

## ----extract_tables, message=FALSE, warning=FALSE-----------------------------
res <- panoramic_extract_contrast(se_meta)
head(res[, c("ct1", "ct2", "radius_um", "beta_diff", "p_diff", "fdr_diff")])
res %>% dplyr::filter(fdr_diff < 0.05)

## ----one_call, message=FALSE, warning=FALSE-----------------------------------
out <- panoramic_analyze(
  spe_list = spe_list,
  design = design,
  cell_type = "cell_type",
  radii_um = radii_um,
  nsim = 20,
  min_cells = 5,
  window = "concave",
  BPPARAM = BiocParallel::SerialParam()
)
names(out)
names(out$tables)

## ----forest_plot, message=FALSE, warning=FALSE--------------------------------
plot_forest(
  se_meta,
  ct1 = "A",
  ct2 = "B",
  radius_um = 25,
  group_col = "group"
)

## ----volcano_plot, message=FALSE, warning=FALSE-------------------------------
plot_volcano(se_meta)

## ----network_plot, message=FALSE, warning=FALSE-------------------------------
net <- create_spatial_network(
  se_meta,
  fdr_threshold = 0.2,
  directed = FALSE,
  leiden_resolution = 1.0
)
plot_spatial_network(
  se_meta,
  fdr_threshold = 0.2,
  directed = FALSE,
  layout = "fr",
  node_size_by = "degree"
)

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