---
title: "Supplementary Code 2 | Post-classification homolog annotation and density-level assignment"
format:
  html:
    toc: true
    toc-depth: 3
    code-fold: false
execute:
  echo: true
  warning: true
  message: true
editor: visual
---

# Purpose

This Quarto document first selects the optimal homolog-classification path from the set of candidate paths generated by the energy-based workflow. The upstream procedure constructs one greedy similarity path from each possible starting nucleus and records the cumulative energy associated with each candidate path. For the chromosome and cumulative-energy variable specified in the analysis settings (by default, `Chr12` and `CumulativeTotalEnergy`), this document extracts the terminal cumulative energy of every candidate path and selects the path with the lowest value. The selected path therefore represents the most spatially coherent candidate path under the original cumulative-energy criterion; no energy calculation, path-construction rule, or homolog-assignment rule is recomputed or modified here.

The records associated with the selected minimum-energy path are then used to annotate chromosome-territory observations with the inferred homolog labels (`Alpha` and `Beta`). Specifically, the script matches the selected-path assignments to the chromosome-territory dataset using the cell identifier and chromosome object name, thereby transferring the inferred homolog identity to each observation. It subsequently assigns within-group kernel-density levels, optionally applies a display-oriented Alpha/Beta label swap, generates an example KDE plot, and saves the annotated dataset.

# Required packages

```{r}
required_packages <- c("rio", "dplyr", "ks", "plotrix", "tcltk")

missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]

if (length(missing_packages) > 0L) {
  stop(
    "Install the following required R packages before running this document: ",
    paste(missing_packages, collapse = ", "),
    call. = FALSE
  )
}
```

# Import input datasets

## Chromosome-territory dataset

```{r}
# Select the Imaris-derived chromosome-territory dataset.
chromosome_file <- tcltk::tk_choose.files(
  caption = "Select the chromosome-territory dataset"
)

if (length(chromosome_file) != 1L || !nzchar(chromosome_file)) {
  stop("Exactly one chromosome-territory dataset must be selected.", call. = FALSE)
}

# Keep the original object name for compatibility with the analysis workflow.
DataSet_Cromosomi <- rio::import(chromosome_file)

rm(chromosome_file)
```

## Energy-classification dataset

```{r}
# Select the RData/RDS/importable object produced by the homolog-classification workflow.
energy_file <- tcltk::tk_choose.files(
  caption = "Select the energy-classification dataset"
)

if (length(energy_file) != 1L || !nzchar(energy_file)) {
  stop("Exactly one energy-classification dataset must be selected.", call. = FALSE)
}

Dati_Energia <- rio::import(energy_file)

rm(energy_file)
```

# Compatibility layer for energy-result objects

```{r}
#' Standardize energy-result element names without changing their content.
#'
#' Accepted input conventions:
#'   1. lista / dataframe
#'   2. paths / data
#'
#' Returned convention:
#'   $lista     : one energy path per starting cell
#'   $dataframe : concatenated observations across all paths
standardize_energy_output <- function(energy_results) {
  if (!is.list(energy_results)) {
    stop("The imported energy result must be a list.", call. = FALSE)
  }

  input_names <- names(energy_results)

  path_element <- if ("lista" %in% input_names) {
    "lista"
  } else if ("paths" %in% input_names) {
    "paths"
  } else {
    NA_character_
  }

  data_element <- if ("dataframe" %in% input_names) {
    "dataframe"
  } else if ("data" %in% input_names) {
    "data"
  } else {
    NA_character_
  }

  if (is.na(path_element) || is.na(data_element)) {
    stop(
      "Energy results must contain either 'lista'/'dataframe' or 'paths'/'data'. ",
      "Available elements: ", paste(input_names, collapse = ", "),
      call. = FALSE
    )
  }

  list(
    lista = energy_results[[path_element]],
    dataframe = energy_results[[data_element]]
  )
}

# From this point onward, the analysis always uses the canonical names.
Dati_Energia <- standardize_energy_output(Dati_Energia)
```

# Path selection and homolog-label annotation

```{r}
#' Select the globally lowest- and highest-energy paths.
#'
#' For each starting cell, the terminal row of the selected chromosome-specific
#' path is extracted. The function then returns the path(s) with the global
#' minimum and maximum value of the requested cumulative-energy variable.
#'
#' This preserves the original path-selection logic.
Path_min_max_energia_AlphaBeta <- function(
  lista_dati,
  oggetto,
  variabile_min_max
) {
  terminal_rows <- lapply(lista_dati[["lista"]], function(path_by_chromosome) {
    path_by_chromosome[[oggetto]] |>
      dplyr::slice_tail(n = 1L) |>
      dplyr::filter(
        .data[[variabile_min_max]] == min(.data[[variabile_min_max]])
      )
  })

  terminal_energy_table <- as.data.frame(do.call(rbind, terminal_rows))

  Path_ottimale <- terminal_energy_table |>
    dplyr::filter(
      .data[[variabile_min_max]] == min(.data[[variabile_min_max]])
    ) |>
    dplyr::select(
      Object_Name,
      Pattern,
      Path,
      dplyr::all_of(variabile_min_max)
    )

  Path_sfavorevole <- terminal_energy_table |>
    dplyr::filter(
      .data[[variabile_min_max]] == max(.data[[variabile_min_max]])
    ) |>
    dplyr::select(
      Object_Name,
      Pattern,
      Path,
      dplyr::all_of(variabile_min_max)
    )

  list(
    migliore = Path_ottimale,
    peggiore = Path_sfavorevole
  )
}

#' Append energy-derived Alpha/Beta labels to chromosome-territory observations.
#'
#' The joining keys are IdentificativoCellula and Object_Name, exactly as in the
#' original workflow. The input label equal to Alpha is annotated as "Alpha";
#' all remaining labels are annotated as "Beta", preserving the original rule.
process_data <- function(lista_dati, dati_chr, Path_ottimo) {
  selected_energy_records <- lista_dati[["dataframe"]] |>
    dplyr::filter(Path == Path_ottimo)

  homolog_lookup <- selected_energy_records |>
    dplyr::select(
      IdentificativoCellula,
      Object_Name,
      Frame,
      Alpha,
      Beta
    )

  annotated_records <- dplyr::left_join(
    dati_chr,
    homolog_lookup,
    by = c("IdentificativoCellula", "Object_Name")
  ) |>
    dplyr::mutate(
      Omologo_Energetico = ifelse(
        Omologo == Alpha,
        "Alpha",
        "Beta"
      )
    )

  annotated_records |>
    dplyr::select(
      IdentificativoCellula,
      Frame,
      Pattern,
      Orientamento,
      Replica,
      Object_Name,
      Omologo,
      Omologo_Energetico,
      X_circ,
      Y_circ,
      Z_circ,
      Raggio_Circ,
      Angolo_Circ,
      X_cart,
      Y_cart,
      Z_cart,
      Raggio,
      Angolo,
      Dist.Chr2Perim.Norm
    )
}
```

# Select the energetically optimal path

```{r}
# Analysis settings. Modify these values only when analyzing another condition.
analysis_pattern <- "F0.M0.0"
path_selection_chromosome <- "Chr12"
energy_variable <- "CumulativeTotalEnergy"

# Identify the path with the lowest terminal cumulative energy.
path_summary <- Path_min_max_energia_AlphaBeta(
  lista_dati = Dati_Energia,
  oggetto = path_selection_chromosome,
  variabile_min_max = energy_variable
)

Path_ottimo <- path_summary[["migliore"]][["Path"]]

# Restrict chromosome-territory records to the selected experimental condition.
dati_chr <- DataSet_Cromosomi |>
  dplyr::filter(Pattern == analysis_pattern)

# Add energy-derived Alpha/Beta labels to the selected condition.
dataset <- process_data(
  lista_dati = Dati_Energia,
  dati_chr = dati_chr,
  Path_ottimo = Path_ottimo
)

rm(path_summary, dati_chr)
```

# Density-level assignment

```{r}
#' Assign density quartiles within one chromosome/homolog subgroup.
#'
#' The procedure is identical to the original repeated code blocks:
#'   1. fit a two-dimensional KDE to X_circ and Y_circ;
#'   2. predict density at each observed coordinate;
#'   3. assign levels 1–4 using the 75th, 50th, and 25th percentiles.
#'
#' Level 1 corresponds to the highest-density quartile and level 4 to the
#' lowest-density quartile.
assign_density_levels <- function(data, object_name, homolog_label) {
  subgroup <- data |>
    dplyr::filter(
      Object_Name %in% c(object_name),
      Omologo_Energetico %in% c(homolog_label)
    )

  kde_coordinates <- subgroup |>
    dplyr::select(X_circ, Y_circ)

  density_model <- ks::kde(kde_coordinates)
  estimated_density <- predict(density_model, x = kde_coordinates)

  subgroup_with_density <- subgroup |>
    dplyr::mutate(density = estimated_density)

  density_quartiles <- quantile(
    subgroup_with_density$density,
    probs = c(0.25, 0.50, 0.75)
  )

  subgroup_with_density |>
    dplyr::mutate(
      Livello = dplyr::case_when(
        density >= density_quartiles[3] ~ 1,
        density >= density_quartiles[2] &
          density < density_quartiles[3] ~ 2,
        density >= density_quartiles[1] &
          density < density_quartiles[2] ~ 3,
        density < density_quartiles[1] ~ 4
      )
    )
}

# Preserve the original output object names used by downstream analyses.
dati_Chr06_Alpha <- assign_density_levels(dataset, "Chr06", "Alpha")
dati_Chr06_Beta  <- assign_density_levels(dataset, "Chr06", "Beta")
dati_Chr12_Alpha <- assign_density_levels(dataset, "Chr12", "Alpha")
dati_Chr12_Beta  <- assign_density_levels(dataset, "Chr12", "Beta")

dataset_livelli <- dplyr::bind_rows(
  dati_Chr06_Alpha,
  dati_Chr06_Beta,
  dati_Chr12_Alpha,
  dati_Chr12_Beta
)
```

# Optional label calibration for display

```{r}
#' Swap energy-derived Alpha/Beta labels for one chromosome.
#'
#' This operation affects only Omologo_Energetico and is retained from the
#' original workflow for display-oriented calibration.
scambia_nomi_livelli <- function(df, object_name) {
  df |>
    dplyr::mutate(
      Omologo_Energetico = dplyr::case_when(
        Object_Name == object_name &
          Omologo_Energetico == "Alpha" ~ "Beta",
        Object_Name == object_name &
          Omologo_Energetico == "Beta" ~ "Alpha",
        TRUE ~ Omologo_Energetico
      )
    )
}

dataset <- dataset_livelli

# Retained original calibration: swap labels for chromosome 6 only.
#dataset <- scambia_nomi_livelli(dataset, "Chr06")

# To swap labels for chromosome 12 instead, use:
dataset <- scambia_nomi_livelli(dataset, "Chr12")
```

# Example KDE plot

```{r}
# Select the chromosome and energy-derived homolog to visualize.
obj_name <- "Chr12"
omolog_plot <- "Beta"

dati_KDE <- dataset |>
  dplyr::filter(
    Object_Name %in% c(obj_name),
    Omologo_Energetico %in% c(omolog_plot)
  ) |>
  dplyr::select(X_circ, Y_circ)

densità <- ks::kde(dati_KDE)

# Select the same colour palettes used in the original workflow.
palette_colori <- if (obj_name == "Chr06") {
  c("white", "#FF9B9B", "#F24C3D", "#DC0000", "#850000", "#850000")
} else if (obj_name == "Chr12") {
  c("white", "#E6FF94", "#9DDE8B", "#40A578", "#006769", "#006769")
} else {
  c("white", "lightgray", "gray", "darkgray", "black", "black")
}

par(
  mar = c(3, 3, 0.5, 0.5),
  mgp = c(1.8, 0.5, 0),
  cex.axis = 1,
  cex.lab = 1.1,
  bty = "n",
  pty = "s"
)

plot(
  densità,
  display = "filled.contour",
  drawlabels = FALSE,
  cont = seq(0, 100, by = 25),
  col = palette_colori,
  xlab = expression("X (a.u.)"),
  ylab = expression("Y (a.u.)"),
  xlim = c(-1, 1),
  ylim = c(-1, 1),
  axes = FALSE,
  asp = 1,
  frame.plot = FALSE,
  main = NULL
)

plotrix::draw.ellipse(
  0,
  0,
  a = 1,
  b = 1,
  nv = 100,
  border = "#003285",
  lwd = 1.5,
  lty = 2
)

axis(side = 1, at = seq(-1, 1, 1), cex.axis = 1, tcl = -0.2, lwd = 1.5)
axis(side = 2, at = seq(-1, 1, 1), cex.axis = 1, tcl = -0.2, lwd = 1.5)

rm(dati_KDE, obj_name)
```

# Save the annotated dataset

```{r}
output_directory <- tcltk::tk_choose.dir(
  caption = "Select the directory in which to save the annotated dataset"
)

if (length(output_directory) != 1L || !nzchar(output_directory)) {
  stop("An output directory must be selected.", call. = FALSE)
}

save(
  dataset,
  file = file.path(output_directory, "02_Demo_Top_Postclassification.RData")
)

rm(output_directory)
```
