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

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

## ----example-data-------------------------------------------------------------
library(CySA)
sce <- CySA_example_sce()
head(S4Vectors::metadata(sce)$SOM_codes)

## ----prep-data----------------------------------------------------------------
prepped <- prepClusterSelectorData(
  sce,
  total_cells_to_sample = 200,
  somCodesName = "SOM_codes"
)
names(prepped)

## ----build-inputs-------------------------------------------------------------
som_codes <- S4Vectors::metadata(sce)$SOM_codes
markers <- S4Vectors::metadata(sce)$map$colsUsed

dend <- stats::as.dendrogram(stats::hclust(stats::dist(som_codes)))

dendTable <- data.frame(
  id = seq_len(nrow(som_codes)),
  label = rownames(som_codes),
  stringsAsFactors = FALSE
)

clusterPatientTable <- table(
  sample_id = sce$sample_id,
  cluster_id = sce$cluster_id
)

somRasterData <- data.frame(
  x = rep(seq_len(5), length.out = nrow(som_codes)),
  y = rep(seq_len(2), each = ceiling(nrow(som_codes) / 2)),
  id = seq_len(nrow(som_codes))
)
for (m in markers) {
  somRasterData[[m]] <- seq_len(nrow(som_codes)) / nrow(som_codes)
}

arr <- array(
  data = seq_len(10 * 10 * length(markers)),
  dim = c(10, 10, length(markers))
)
somRasterObj <- raster::brick(arr)
names(somRasterObj) <- markers

## ----cluster-selector---------------------------------------------------------
app <- clusterSelector(
  sce = prepped$sce,
  sce_subsampled = prepped$sce_subsampled,
  dList = prepped$dList,
  dend = dend,
  dendTable = dendTable,
  clusterPatientTable = clusterPatientTable,
  somRasterData = somRasterData,
  somRasterObj = somRasterObj
)

## ----launch-app, eval = FALSE-------------------------------------------------
# if (interactive()) {
#   shiny::runApp(app)
# }

## ----fig-gating, echo = FALSE, out.width = "49%", fig.show = "hold"-----------
knitr::include_graphics(c(
  "../man/figures/README-2d-plot.png",
  "../man/figures/README-dendrogram.png"
))

## ----fig-dimred, echo = FALSE, out.width = "49%", fig.show = "hold"-----------
knitr::include_graphics(c(
  "../man/figures/README-som-2d-plots.png",
  "../man/figures/README-dimreduction.png"
))

## ----flowsom-prep, eval = FALSE-----------------------------------------------
# library(flowCore)
# library(flowWorkspace)
# library(FlowSOM)
# 
# # Load example FCS file (shipped with FlowSOM)
# fileName <- system.file("extdata", "68983.fcs", package = "FlowSOM")
# ff <- read.FCS(fileName)
# 
# # Compensation (for flow cytometry data)
# comp <- keyword(ff)[["SPILL"]]
# ff <- compensate(ff, comp)
# 
# # Transformation
# # For flow cytometry: logicle transformation
# transformList <- estimateLogicle(ff, channels = colnames(comp))
# ff <- transform(ff, transformList)
# 
# # For CyTOF data, use arcsinh transformation instead:
# # ff <- transform(ff, arcsinhTransform())

## ----flowsom-run, eval = FALSE------------------------------------------------
# set.seed(42)
# 
# # Run FlowSOM on the preprocessed flowFrame
# fSOM <- FlowSOM(ff,
#   # Input options:
#   compensate = FALSE,    # already compensated above
#   transform = FALSE,     # already transformed above
#   scale = FALSE,
#   # SOM options:
#   colsToUse = c(9, 12, 14:18),  # select relevant channels
#   xdim = 7, ydim = 7,           # 7x7 SOM grid = 49 nodes
#   # Metaclustering options:
#   nClus = 10                    # number of metaclusters
# )
# 

## ----flowsom-clusters, eval = FALSE-------------------------------------------
# head(GetClusters(fSOM))      # SOM node ID for each cell
# head(GetMetaclusters(fSOM))  # Metacluster ID for each cell

## ----flowsom-to-sce, eval = FALSE---------------------------------------------
# library(SingleCellExperiment)
# library(S4Vectors)
# 
# # Extract SOM codes (one row per SOM node, one column per marker)
# som_codes <- fSOM$map$codes
# markers <- fSOM$map$colsUsed  # Get the marker names used in the SOM
# 
# # Extract cell-level assignments
# cell_clusters <- GetClusters(fSOM)
# cell_metaclusters <- GetMetaclusters(fSOM)
# 
# # Build colData from the flowFrame
# # CRITICAL: cluster_id must be a factor with levels for ALL SOM nodes (1 to nrow(som_codes)),
# # even if some nodes are empty. Otherwise, prepClusterSelectorData() will fail when it
# # tries to set levels based on nrow(SOM_codes), and plot interactions will break because
# # node IDs won't map correctly.
# coldata <- data.frame(
#   sample_id = rep("sample1", nrow(ff@exprs)),
#   cluster_id = factor(cell_clusters, levels = seq_len(nrow(som_codes))),
#   metacluster_id = factor(cell_metaclusters)
# )
# 
# # Create the expression matrix with ONLY the markers used in the SOM
# # CRITICAL: The rownames of the expression matrix MUST match the column names of som_codes.
# # A common mistake is to use all flowFrame channels (colnames(ff)), which includes FSC, SSC,
# # Time, and other parameters not used in the SOM. When CySA's plotSOMScatter() tries to
# # index som_codes with marker names that don't exist, you get errors like:
# #   "subscript out of bounds" or "dim(X) must have a positive length"
# # The fix: subset to only the markers that were actually used in the FlowSOM analysis.
# exprs_mat <- t(ff@exprs[, markers, drop = FALSE])
# rownames(exprs_mat) <- markers
# 
# # Compute SOM_stats (required by CySA)
# # CRITICAL: SOM_stats must have one row per SOM node, with an 'id' column.
# # A common mistake is to compute statistics per marker instead of per node.
# # The 'id' values must be 1:nrow(som_codes) to match the SOM node indices.
# som_stats <- data.frame(
#   id = seq_len(nrow(som_codes)),
#   n = tabulate(cell_clusters, nbins = nrow(som_codes)),
#   mean = som_codes[, 1],  # Use first marker code as proxy for visualization
#   median = som_codes[, 1],
#   rdQu = som_codes[, 1],
#   max = som_codes[, 1],
#   stringsAsFactors = FALSE
# )
# rownames(som_stats) <- rownames(som_codes)
# 
# # Create the SingleCellExperiment object
# # CRITICAL: experiment_info must contain at least one NUMERIC column (besides sample_id).
# # The Stats panel tries to select numeric columns from experiment_info to offer as
# # normalization options. If there are no numeric columns, apply() fails with:
# #   "dim(X) must have a positive length"
# # The fix: include at least one numeric column like total_cells or sample_nr.
# experiment_info <- data.frame(
#   sample_id = unique(coldata$sample_id),
#   total_cells = nrow(ff@exprs),  # numeric column required by stats panel
#   stringsAsFactors = FALSE
# )
# 
# sce <- SingleCellExperiment(
#   assays = list(exprs = exprs_mat),
#   colData = DataFrame(coldata),
#   metadata = list(
#     SOM_codes = som_codes,
#     SOM_stats = som_stats,
#     map = list(colsUsed = markers),
#     experiment_info = experiment_info
#   )
# )
# 
# # Now pass to CySA
# # CRITICAL: prepClusterSelectorData() requires at least 12 markers to build the default
# # dList (marker pairs for 2D plots). FlowSOM often uses fewer markers (e.g., 7 in this
# # example). If you get the error:
# #   "sce must have at least 12 row names to build default dList"
# # The fix: provide a custom dList with explicit marker pairs.
# dList <- list(
#   d1 = c(markers[1], markers[2]),
#   d2 = c(markers[3], markers[4]),
#   d3 = c(markers[5], markers[6]),
#   d4 = c(markers[1], markers[3]),
#   d5 = c(markers[2], markers[4]),
#   d6 = c(markers[5], markers[7])
# )
# prepped <- prepClusterSelectorData(
#   sce,
#   total_cells_to_sample = 500,
#   dList = dList
# )

## ----flowsom-app, eval = FALSE------------------------------------------------
# # Extract SOM codes and build dendrogram
# som_codes <- fSOM$map$codes
# 
# # CRITICAL: FlowSOM may not set rownames on the codes matrix.
# # Without rownames, the dendTable will have mismatched row counts, causing:
# #   "arguments imply differing number of rows"
# # The fix: assign node IDs as rownames if they're missing.
# if (is.null(rownames(som_codes))) {
#   rownames(som_codes) <- seq_len(nrow(som_codes))
# }
# 
# # Build hierarchical clustering dendrogram from SOM codes
# # This is used for the dendrogram view in CySA
# dend <- stats::as.dendrogram(stats::hclust(stats::dist(som_codes)))
# 
# # Build the dendrogram navigation table
# # CRITICAL: The 'id' column must match the SOM node indices (1:nrow(som_codes))
# # and 'label' must match the rownames of som_codes.
# dendTable <- data.frame(
#   id = seq_len(nrow(som_codes)),
#   label = rownames(som_codes),
#   stringsAsFactors = FALSE
# )
# 
# # Build cluster-by-sample table (SOM nodes x samples)
# # This is used for abundance comparisons across samples
# clusterPatientTable <- table(
#   sample_id = sce$sample_id,
#   cluster_id = sce$cluster_id
# )
# 
# # Build SOM raster data for heatmap visualization
# # CRITICAL: The SOM grid layout must match the FlowSOM dimensions.
# # FlowSOM uses a rectangular grid (xdim × ydim). The raster data frame
# # must have columns: x, y, id, and one column per marker with the SOM codes.
# # A common mistake is to use the wrong grid dimensions, which causes:
# #   "subscript out of bounds" or misaligned heatmaps
# xdim <- fSOM$map$xdim
# ydim <- fSOM$map$ydim
# 
# # Create grid coordinates for each SOM node
# # Nodes are arranged row-by-row in the FlowSOM grid
# somRasterData <- data.frame(
#   x = rep(seq_len(xdim), length.out = nrow(som_codes)),
#   y = rep(seq_len(ydim), each = ceiling(nrow(som_codes) / ydim))[seq_len(nrow(som_codes))],
#   id = seq_len(nrow(som_codes))
# )
# 
# # Add marker expression values for each SOM node
# # These are used to color the SOM heatmap tiles
# for (m in colnames(som_codes)) {
#   somRasterData[[m]] <- som_codes[, m]
# }
# 
# # Create the CySA app
# # CRITICAL: All arguments must be provided (no NULLs allowed for required inputs).
# # The most common errors at this stage are:
# #   - Missing experiment_info numeric column → "dim(X) must have a positive length"
# #   - Mismatched marker names between exprs and som_codes → "subscript out of bounds"
# #   - cluster_id factor missing levels → "undefined factor levels"
# app <- clusterSelector(
#   sce = prepped$sce,
#   sce_subsampled = prepped$sce_subsampled,
#   dList = prepped$dList,
#   dend = dend,
#   dendTable = dendTable,
#   clusterPatientTable = clusterPatientTable,
#   somRasterData = somRasterData,
#   somRasterObj = NULL  # not needed when somRasterData is provided
# )
# 
# # Launch interactively
# if (interactive()) {
#   shiny::runApp(app)
# }

## ----fig-stats, echo = FALSE, out.width = "100%"------------------------------
knitr::include_graphics("../man/figures/README-stats-panel.png")

## ----example with stats-------------------------------------------------------
library(CySA)
library(SingleCellExperiment)
library(S4Vectors)
set.seed(42)

# ── dimensions ────────────────────────────────────────────────────────────────
n_markers <- 12
n_som_nodes <- 50
n_samples <- 6 # 3 control + 3 treated
n_cells <- 300 # per sample

marker_names <- paste0("marker", seq_len(n_markers))
sample_ids <- paste0("S", seq_len(n_samples))
conditions <- factor(
  c(rep("control", 3), rep("treated", 3)),
  levels = c("control", "treated")
)

# ── experiment_info ───────────────────────────────────────────────────────────
# total_cells: cells acquired by the cytometer (used as relativeTo denominator)
experiment_info <- data.frame(
  sample_id = sample_ids,
  condition = conditions,
  total_cells = c(8000L, 9200L, 7800L, 10500L, 11000L, 9800L),
  stringsAsFactors = FALSE
)

# ── SOM codes ─────────────────────────────────────────────────────────────────
# Nodes  1-20 : low expression  → "resting"  phenotype
# Nodes 21-30 : intermediate
# Nodes 31-50 : high expression → "activated" phenotype
som_codes <- matrix(0,
  nrow = n_som_nodes,
  ncol = n_markers,
  dimnames = list(
    paste0("node", seq_len(n_som_nodes)),
    marker_names
  )
)
for (node in seq_len(n_som_nodes)) {
  base <- node / n_som_nodes # 0.02 … 1.00
  som_codes[node, ] <- pmax(0, base + rnorm(n_markers, sd = 0.02))
}

# ── cell assignment ───────────────────────────────────────────────────────────
# Control: 70 % in nodes  1-20  (resting)
# Treated: 70 % in nodes 31-50  (activated)
assign_clusters <- function(condition, n) {
  if (condition == "control") {
    c(
      sample(1:20, round(n * 0.70), replace = TRUE),
      sample(seq_len(n_som_nodes), n - round(n * 0.70), replace = TRUE)
    )
  } else {
    c(
      sample(31:50, round(n * 0.70), replace = TRUE),
      sample(seq_len(n_som_nodes), n - round(n * 0.70), replace = TRUE)
    )
  }
}

coldata_list <- mapply(function(sid, cond) {
  data.frame(
    sample_id = sid,
    cluster_id = assign_clusters(as.character(cond), n_cells),
    stringsAsFactors = FALSE
  )
}, sample_ids, as.character(conditions), SIMPLIFY = FALSE)

coldata_df <- do.call(rbind, coldata_list)
n_total <- nrow(coldata_df) # 1800

# ── assay matrix — each cell ≈ its node's SOM code + noise ───────────────────
exprs_mat <- vapply(seq_len(n_total), function(i) {
  node <- coldata_df$cluster_id[i]
  pmax(0, som_codes[node, ] + rnorm(n_markers, sd = 0.05))
}, numeric(n_markers))
dimnames(exprs_mat) <- list(marker_names, paste0("cell", seq_len(n_total)))

# ── SOM_stats ─────────────────────────────────────────────────────────────────
node_counts <- tabulate(coldata_df$cluster_id, nbins = n_som_nodes)

# summarise marker1 expression per node for hover-text columns
m1 <- exprs_mat["marker1", ]
node_factor <- factor(coldata_df$cluster_id, levels = seq_len(n_som_nodes))

som_stats <- data.frame(
  id = seq_len(n_som_nodes),
  n = node_counts,
  mean = as.numeric(tapply(m1, node_factor, mean)),
  median = as.numeric(tapply(m1, node_factor, median)),
  rdQu = as.numeric(tapply(m1, node_factor, quantile, probs = 0.75)),
  max = as.numeric(tapply(m1, node_factor, max)),
  stringsAsFactors = FALSE
)
# nodes with zero cells get NA from tapply; replace with 0
som_stats[is.na(som_stats)] <- 0

# ── assemble SCE ──────────────────────────────────────────────────────────────
sce_stats <- SingleCellExperiment(
  assays = list(exprs = exprs_mat),
  colData = DataFrame(
    sample_id  = coldata_df$sample_id,
    cluster_id = factor(coldata_df$cluster_id, levels = seq_len(n_som_nodes))
  ),
  metadata = list(
    SOM_codes       = som_codes,
    SOM_stats       = som_stats,
    map             = list(colsUsed = marker_names),
    experiment_info = experiment_info
  )
)
rownames(sce_stats) <- marker_names

# ── clusterPatientTable ───────────────────────────────────────────────────────
clusterPatientTable <- table(
  sample_id  = sce_stats$sample_id,
  cluster_id = sce_stats$cluster_id
)

# ── app inputs ────────────────────────────────────────────────────────────────
prepped <- prepClusterSelectorData(
  sce_stats,
  total_cells_to_sample = 600,
  somCodesName = "SOM_codes"
)

som_codes_mat <- metadata(sce_stats)$SOM_codes
dend <- stats::as.dendrogram(stats::hclust(stats::dist(som_codes_mat)))

dendTable <- data.frame(
  id = seq_len(n_som_nodes),
  label = rownames(som_codes_mat),
  stringsAsFactors = FALSE
)

# SOM raster grid: 10 × 5 layout for 50 nodes
somRasterData <- data.frame(
  x  = rep(seq_len(10), times = 5),
  y  = rep(seq_len(5), each = 10),
  id = seq_len(n_som_nodes)
)
for (m in marker_names) {
  somRasterData[[m]] <- som_codes_mat[, m]
}

app <- clusterSelector(
  sce                  = prepped$sce,
  sce_subsampled       = prepped$sce_subsampled,
  dList                = prepped$dList,
  dend                 = dend,
  dendTable            = dendTable,
  clusterPatientTable  = clusterPatientTable,
  somRasterData        = somRasterData,
  somRasterObj         = NULL
)

## ----fig-phenotype, echo = FALSE, out.width = "100%"--------------------------
knitr::include_graphics(c(
  "../man/figures/README-marker-pies.png",
  "../man/figures/README-som-pairs-grid.png",
  "../man/figures/README-som-heatmaps.png"
))

## ----launch-app-2, eval = FALSE-----------------------------------------------
# if (interactive()) {
#   shiny::runApp(app)
# }

## ----som-scatter--------------------------------------------------------------
plotSOMScatter(sce, chs = c("marker1", "marker2"))

## ----scatter-bj---------------------------------------------------------------
plotCytoScatter(sce, chs = c("marker1", "marker2"))

## ----session-info-------------------------------------------------------------
sessionInfo()